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ABSTRACT 


A  classic  problem  in  the  study  of  the  physics  of  sound 
in  fluids  is  that  of  the  finite  amplitude  plane  wave  propaga- 
ting in  a  viscous,  unbounded  medium.   Though  the  solution  to 
this  problem  is  well  known,  for  a  given  range  of  parameters, 
it  would  be  desirable  to  develop  techniques  for  solution   of 
the  problem  which  are  not  similarly  limited.   Through  the  ap- 
plication of  the  technique  of  parametric  differentiation, 
this  goal  is  realizable.   This  extension  method  transforms 
the  governing  non-linear  differential  equation  to  a  linear 
equation  in  what  may  be  termed  parameter  space;  the  equation 
is  solved  and  a  quadrature  recovers  the  dependent  variable 
solution.  Parametric  differentiation  is  conceptually  straight- 
forward and  has  been  applied  in  the  past  to  a  wide  variety  of 
non-linear  equations.   The  method  has  b'  n  applied  to  the  e- 
quation  which  describes  finite  amplitude  plane  wave  propaga- 
tion in  a  viscous  fluid  in  order  to  compare  the  results  to  the 
predictions  of  a  known  perturbation  solution;  the  ultimate  ob- 
jective being  to  utilize  the  technique  for  exact  solution  of 
the  corresponding  spherical  and  cylindrical  wave  problems. 
The  first  step  in  this  process  has  been  achieved;  that  is, 
numerical  solutions  for  the  plane  wave  case  have  been  gener- 
ated for  a  range  of  values  of  the  various  viscous  and  non- 
linear parameters  which  are  consistent  with  results  obtained 
analytically.   Further  it  is  shown  that  solutions  generated 
through  the  application  of  parametric  differentiation  may  in 
fact  have  greater  validity  for  certain  ranges  of  viscous  and 
non-linear  parameters. 

Thesis  Supervisor:   Dr.  Wesley  L.  Harris,  Sr. 

Title:  Associate  Professor  of  Ocean  Engineering 
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NOTATION 


A       the  coefficient  of  the  first-order  term  in  an  assumed 
liquid  pressure-density  relation 

A,,  A9  ,  A.,,  A.  -  the  variable  coefficients  of  g   ,  g   ,  g  , 

and  g     respectively 

a,  c     arbitrary  coefficients  in  a  generalized  boundary  con- 
dition 

2 
a        small  signal  attenuation  coefficient  =  vbk  /2C   for 

liquids,  =  v[Y  +  (y-1) /Pr] k2/2C   for  gases 

B       the  coefficient  of  the  second-order  term  in  an  assumed 
liquid  pressure-density  relation 


a  viscosity  number  for  liquids  which  represents  shear 
viscosity  and  a  phenomonological  bulk  viscosity 


3       the  parameter  of  non-linearity  =  1  +  B/2A  for  liquids, 
=  (y+l)/2  for  gases 


COO,  C01 ,  etc.  -  coefficients  of  the  points  of  the  finite-dif- 
ference cube 


C       small  signal  sound  speed 


C       specific  heat  at  constant  pressure 


C       specific  heat  at  constant  volume 


interval  or  change  in  the  variable  which  it  precedes 
acoustic  Mach  number 


e       Neumann  factor 
n 
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n,  n  '    coefficients  of  shear  and  bulk  viscosity  respectively 


the  variable  inhomogeneous  term  of  the  parameter  space 
equation 


parameter  space  variable  =  9£/8ip 


9^'  <3„„>    9^4-'  Svv'  etc-  "  notation  for  dg/dx,    8  2g/8x2  ,  8  3g/8x2  8t , 

82^/8x2,  etc. 


X    ^XX'   ^XXt     XX  2r/~> 


g.,  E.   values  of  these  variables  for  the  ith  value  of  \b ,    \\) . 

g.  .     indicial  notation  for  the  two-dimensional  solution  grid 

T       an  indicator  of  the  strength  of  the  non-linearity  rela- 
tive to  dissipation  =  $ek/a 


Y  ratio  of  specific  heats  =  C  /C 
'  r  p'  v 

I  nth  order  modified  Bessel  function  of  the  first  kind 

n  — 

J  nth  order  Bessel  function  of  the  first  kind 

n  — 

k  acoustic  wave  number 


k        3£k  product 

L[  ]     a  general  linear  operator 


wavelength  of  acoustic  signal,  also  used  as  an  arbitrary 
real  number  in  VonNeumann  stability  test 


m,  n     indicial  notation  as  in  g    ,  corresponds  to  g . 


N[  ]     a  general  non-linear  operator 
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v       kinematic  viscosity  =  n/p 

5       particle  displacement  amplitude  =  U  /oj 

£,  £     particle  displacement;  subscript  indicates  the  base 
solution,  i=0 


£f      an  expression  for  a  calculated  particle  displacement 
for  output 


co       angular  frequency 


P,  Q    constants  in  the  fluid  equation  of  state,  for  gases 


o 


2 
Mo^c 

and  y  is  determined  empirically 


P  =  PQ,  Q  =  0;  for  liquids,  P  =  P0CqVy,  Q  =  P  -  PQ , 


instantaneous  pressure 


p       ambient  pressure 
*o  r 


Pr    •   Prandtl  number 


<J>,  $n        an  arbitrary  dependent  variable,  subscript  indicates 
the  base  solution 


instantaneous  density 


p       ambient  density 

o  J 


s       condensation  =  p-p  /p 

o   o 


a       a  nondimensional  distance  =  $ekx 
T       shock  thickness 
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t,  t    temporal  coordinate,  characteristic  time 


an  arbitrary  complex  number  in  VonNeumann  stability 
test 


B. 

1 


boundary  values  of  <J> 


u,  u    particle  velocity,  subscript  indicates  base  solution 
U       particle  velocity  amplitude 
X       discontinuity  distance  =  1/Bek 


X,  X    a  generalized  spatial  vector,  the  subscript  indicates 


the  boundary 


x,  x    spatial  coordinate,  characteristic  length 


Y       viscosity  number 


i> ,    \\)  parameter  of  interest  =  vb/C  X,    subscript  indicates 

base  value 


V       u/U 

'    o 


a  timelike  variable  in  the  Burgers  equation  solution 


a  non-physical  variable  used  to  transform  the  Burgers 
equation  to  a  heat  equation 
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I.    INTRODUCTION 

The  propagation  of  acoustic  waves  in  fluids  and  certain 
crystalline  substances  has  been  the  subject  of  long  and  ex- 
haustive study.   Due  to  the  inherent  simplifications  which  a- 
rise  in  the  examination  of  plane  wave  propagation,  this  prob- 
lem has  received  the  closest  scrutiny.   In  the  general  linear 
theory  of  sound  propagation,  the  assumption  must  be  made  that 
particle  velocity  amplitudes  are  infinitesimal,  allowing  the 
reduction  of  the  problem  to  the  classic  wave  equation.   How- 
ever, in  many  of  today's  acoustic  systems  this  assumption  is 
no  longer  valid,  and  solution  of  the  fully  non-linear  wave 
propagation  problem  is  required  for  analysis  and  prediction 
for  practical  systems.   It  has  long  been  recognized  that  sound 
waves  whose  particle  velocity  amplitudes  are  not  infinitesimal 
have  phase  velocities  whose  magnitudes  change  with  the  local 
particle  velocity,  a  result  which  is  not  predicted  by  the  lin- 
ear theory.   This  phenomenom  received  some  attention  from  the 
philosophers  of  the  nineteenth  century,  and  with  the  exception 
of  the  efforts  of  two  researchers  of  the  pre-war  decade  was 
not  treated  until  the  last  two  decades  at  which  time  signifi- 
cant interest  was  generated  because  of  the  increasing  sophis- 
tication and  power  levels  of  acoustic  systems. 

The  solutions  of  the  1930*3  form  the  basis  of  this  re- 
search.  Fubini  [1]  derived  an  implicit  solution  to  the  plane- 
wave  problem  in  an  inviscid  fluid  and  reduced  it  to  an  explicit 


-10- 


solution  for  low  Mach  number.   Fay  [2] ,  the  other  researcher 
of  the  30 's,  solved  the  same  equation  with  viscous  effects  in- 
cluded.  Blackstock  [3]  noted  that  though  both  solutions  were 
entirely  correct,  they  were  restricted  in  their  spatial  re- 
gions of  applicability.   Applying  weak-shock  theory,  he  demon- 
strated the  manner  in  which  the  two  solutions  are  related  and 
developed  a  single  function  which  describes  the  solution  for 
all  space  in  the  inviscid  case.   Keck  and  Beyer  [4]  performed 
a  perturbation  analysis  which,  for  a  specific  range  of  para- 
meters, described  propagation  in  the  viscous   case.   Black- 
stock  [5] ,  using  approximations  which  he  demonstrated  [6]  to 
be  equivalent  to  those  used  by  Keck  and  Beyer,  was  able  to  re- 
duce the  governing  differential  equation  for  plane  wave  propa- 
gation in  a  thermoviscous  fluid  to  the  Burgers   equation  for  a 
boundary  value  problem,  which  yielded  an  analytic  solution  for 
all  space  and  time.   Propagation  in  relaxing  fluids  has  re- 
ceived moderate  attention  [7,8],  but  has  generally  been  a- 
voided.   Blackstock  [9]  has  derived  analogous  inviscid  solu- 
tions for  cylindrical  and  spherical  waves,  but  as  fate  would 
have  it,  an  analytic  solution  for  cylindrical  and  spherical 
waves  in  a  viscous  fluid  has  not  yet  been  obtained  due  to  the 
fact  that  there  exists  no  apparent  simple  transform,  such  as 
that  utilized  in  the  plane  wave  case,  to  obtain  the  Burgers 
equation. 

There  are  various  parameters  which  describe  the  proper- 
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ties  of  the  medium  and  two  which  relate  to  the  source  excita- 
tion, and  taken  together  they  give  one  an  indication  of  the 
relative  strengths  of  non-linear  effects  and  viscous  (dissipa- 
tive)  effects.   The  parameters  of  the  medium  are  v,  the  kine- 
matic viscosity;  T,  the  viscosity  number;  y ,    the  ratio  of  spe- 
cific heats;  Pr,  the  Prandtl  number;  $,  the  parameter  of  non- 
linearity;  and  C  ,  the  small  signal  sound  speed.  The  expression 

Y-l 

(^  +  -*=■ — )  is  a  collection  of  thermal  and  viscous  coefficients 
Pr 

which  together  with  kinematic  viscosity  relate  the  thermovis- 
cous  nature  of  gases.   In  liquids,  such  an  expression  is  not 
readily  determined  and  is  generally  replaced  by  a  single  term, 
b.    $  is  equal  to  (y  +  l)/2  for  perfect  gases  and  to  1  +  B/2A 
for  fluids  of  arbitrary  equation  of  state  where  A  is  the  coef- 
ficient of  the  first  order  term  and  B/2  is  the  coefficient  of 
the  second  order  term  in  an  assumed  pressure-density  relation. 
Analogous  parameters  for  dielectric  crystals  allow  them  to  be 
treated  in  a  manner  similar  to  that  for  thermoviscous  fluids 
[10] .    The  two  parameters  related  to  the  source  excitation 
are  U  ,  the  peak  particle  velocity,  and  to,  the  angular  fre- 
quency of  the  source.   Much  of  the  previous  effort  in  this 
field  has  made  use  of  a  parameter  Y   =    3ek/a  to  express  the 
relative  strength  of  non-linearity  and  viscosity,  where  e  = 

U  /C  ,  the  acoustic  Mach  number;  k  =  to/C  ,  the  wave  number; 

2 

and  a  =  vbk  /2C  ,  the  small  signal  attenuation  coefficient. 

T  ~  1  may  be  considered  as  the  borderline  for  the  inception  of 
important  non-linear  effects  [11] . 

Another  important  parameter  is  the  discontinuity  distance 
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X  =  (3ek)   ;  at  this  point,  in  inviscid  propagation,  the  slope 
of  the  particle  velocity  waveform  becomes  negatively  infinite 
due  to  non-linear  generation  of  harmonics,  that  is,  X  indi- 
cates the  point  at  which  a  shock  first  forms.   Shocks  were 
first  predicted  by  the  inviscid  theory,  and  the  discontinuity 
distance  which  has  become  a  reference  point  in  most  non-linear 
theory,  was  predicted  by  the  Fubini  solution,  which,  in  fact, 
is  not  valid  beyond  the  discontinuity  distance.   Shocks  form 
because  the  finite  amplitude  of  the  particle  velocity  causes 
the  compressive  portion  of  the  waveform  to  'catch  up1  to  the 
rarefactive  portion  of  the  waveform;  because  a  multi-valued 
waveform  is  a  physical  impossibility,  a  shock  must  form  unless 
dissipation  reduces  the  wave  amplitude  sufficiently.   Figures 
la  and  lb  show  graphically  such  a  waveform.   The  utility  of  V 
has  been  extended  by  Fenlon  [12]  to  indicate  whether  and  where 
a  shock  will  form  in  the  viscous  case;  for  plane  waves  r  >  4.5 
indicates  that  shock  formation  will  occur,  that  is,  non-linear 
effects  will  be  sufficiently  strong,  so  that  even  with  dissi- 
pation, shocks  will  form.   Shock  formation  marks  the  end  of 
the  first  of  three  identifiable  zones  of  propagation  of  the 
finite  amplitude  wave;  this  is  a  region  of  strong  non-linear 
effects.   The  second  region  of  propagation  is  that  in  which 
non-linear  effects  initially  dominate,  but  gradually  taper  off 
as  absorption  reduces  the  magnitude  of  the  fundamental  and 
even  more  rapidly  the  magnitudes  of  the  non-linearly  generated 
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harmonics.   The  limit  of  this  region  is  reached  when  the  rate 
of  decay  of  the  fundamental  due  to  small  signal  attenuation  is 
matched  by  the  rate  of  decay  due  to  the  generation  of  harmon- 
ics, that  is,  when  the  particle  velocity  amplitudes  are  once 
again  infinitesimal. 

The  aforementioned  solutions  may  be  discussed  with  res- 
pect to  their  relations  to  the  parameters  of  the  problem.  The 
Fubini  solution  predicts  shock  formation  at  the  discontinuity 
distance,  but  the  solution  is  no  longer  convergent  beyond  the 
discontinuity  distance.   Fay  produced  a  solution  for  the  most 
stable  waveform,  a  sawtooth,  to  which  the  finite  amplitude 
wave  generates  beyond  about  3  discontinuity  distances.   The 
Fay  solution  contained  T  as  a  parameter,  but  did  not  reduce  to 
Fubini  for  V   ■>  °°;  Blackstock's  bridging  function  accomplished 
this.   Keck  and  Beyer's  perturbation  solution  is  valid  for 
very  weak  non-linear  effects  (relative  to  dissipative  effects). 
Because  of  the  limited  number  of  terms  generated,  it  is  valid 
for  r  =  0(1)  or  less.   Finally,  Blackstock's  Burgers   equa- 
tion solution  is  valid  for  all  r.   Unfortunately  the  series 
which  represents  the  steady  state  solution  is  very  slowly  con- 
vergent, limiting  the  usefulness  of  the  solution.   Finally,  it 
should  be  noted  that  an  inherent  restriction  on  the  size  of 
the  Mach  number  is  assumed  implicitly  or  explicitly  in  all 

4 

cases,  with  a  maximum  value  of  e  =  .1.   Such  a  value  would  be 
stretching  the  applicability  of  certain  approximations  which 
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have  been  made  in  each  of  the  above  cases. 

In  order  to  dispense  with  the  aforementioned  problems 
with  the  two  viscous  solutions,  and  specifically,  to  essenti- 
ally eliminate  the  approximations  which  were  made  therein,  a- 
nother  technique  must  be  used  which  returns  to  the  basic  par- 
tial differential  equation  and  solves  it  without  approximation. 
Parametric  differentiation  is  such  a  technique.   Given  a  non- 
linear differential  equation  in  which  there  exists  a  parameter, 
one  need  only  assume  that  the  solution  is  continuous  in  that 
parameter  to  a  limiting  value  for  which  a  known  solution  ex- 
ists.  This  technique  has  been  applied  in  the  past  by  Rappert 
and  Landahl  [13]  to  non-linear  flow  problems  and  by  Harris  [14] 
to  aerodynamic  sound  produced  by  a  rotating  cylinder  in  a  vis- 
cous medium.   VanDyke  [15]  has  linked  parametric  differentia- 
tion to  work  he  has  performed  on  the  extension  of  perturbation 
series  solutions.   The  application  of  parametric  differentia- 
tion in  this  case  to  a  problem  whose  approximate  solution  is 
known  will  (1)  shed  some  additional  light  on  the  physics  of 
the  problem,  and  (2)  demonstrate  the  applicability  and  utility 
of  a  new  technique  for  use  in  studies  of  the  propagation  of 
finite  amplitude  cylindrical  and  spherical  waves  in  a  viscous 
medium. 
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II.   THE  EQUATION  OF  MOTION 

A  plane  wave  propagating  isentropically  in  a  homogeneous 
viscous  fluid  may  be  completely  described  by  the  continuity 
equation,  the  Navier-Stokes  equation  for  irrotational  flow  and 
an  equation  of  state.   The  appropriate  equations  are: 


1)   Continuity   p  =  (1  +  3£/3x)p  (2.1) 


2)  Navier-Stokes   p^  ^/9t2  =  -dp/dx   +  (|n  +  n,)834/3x28t 

(2.2) 

3)  State   p  =  P(p/pQ)Y  -  Q   [4]  (2.3) 


where  p  is  the  density,  £  is  the  particle  displacement,  p  is 

the  pressure,  n  is  the  shear  viscosity  coefficient,  n '  is  the 

bulk  viscosity  coefficient,  and  P  and  Q  are  constants  dependent 

on  the  fluid  considered.   The  subscripted  parameters  denote 

rest  values  of  the  respective  quantities.   For  ideal  gases, 

P  =  P  /  Q  =  0,  and  y    is  the  ratio  of  specific  heats,  C  /C  . 

For  liquids,  P  =  p  C  2/y,  where  C   is  the  small  signal  sound 
^     '      o  o  o 

speed,  Q  =  P  -  p  ,  and  y   must  be  determined  empirically  with 
the  assumed  equation  of  state  above.   The  equation  of  state 
for  liquids  is  often  written  in  the  form: 

p  =  p   +  A(£=£a)  +  B/2(£^o.)2  +  o((-^-^)  3)  .  .  .         (2.4) 
°       Po  Po        •  Po   ' 

with  p~Pq  =  s   the  condensation.   For  small  values  of  s,  this 
Po 
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equivalent  to  the  equation  of  state  shown  (2.3).   This  can  be 
shown  by  writing 


P  =  P(l  +  P-P^x)^  _  q  (2.5) 

Po 


and  performing  a  binomial  expansion  for  |- — *-©-|  <  1,  the  result 

Po 

being 


p  =  p  +  yP(PZPa)  +  Y(Y  1)  p(P_Pq)2  + _  Q         (2.6) 

Po        ^       Po 


Then,  with  P  -  Q  =  p  ,  A  is  equivalent  to  yP,  and  B  to  y(y-l)P 
and  the  ratio  B/A  =  y-1.  When  one  uses  this  form  the  implicit 
assumption  is  made  that  terms  of  Ok    Q)  ]  and  greater  are  not 
important.   For  any  fluid  of  low  compressibility  this  is 
clearly  valid  for  small  acoustic  Mach  number.   (Note  that  p/pQ 

Returning  to  Equations  2.1-3,  p  and  p  may  be  eliminated 
to  yield  a  single  equation  in  the  particle  displacement.   From 
Equation  2 . 3 


|£=  EX  (Jl/"1  |£  (2.7) 

3x    p0   p0      dx 

Substituting  Equation  2.1  and  its  x  derivative,  ■—   =  -p  -s — |- 

3   ^  dx  "o9xz 

gr  _2 

(1  +  7T-2-)    ,  into  Equation  2.7,  one  obtains, 
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|E  -  -prd  +  |i,-lrl|l|  (2.8) 


Substituting  in  the  Equation  2.2,  rearranging  terms,  and  in- 
serting P  =  p  C  2/y,  one  obtains  the  following  equation: 


<i  +  i>Y+i#--a&-:0 


4    n  ' 
with  v  =  n/p  /  b  =  3-  +  1L- .   At  this  point  one  may  leave  the 

equation  as  it  stands  and  apply  parametric  differentiation 
directly  or  make  one  additional  simplification  which  resorts 
to  the  assumption  that  the  displacements  though  finite  are 
still  not  large  enough  to  require  the  retention  of  cubic  terms 
in  particle  displacement.  Since  previous  effort  has  been  dir- 
ected along  the  second  path  and  since  this  previous  work  is 
the  point  for  comparison  to  determine  the  validity  of  the  ap- 
proach,  this  additional  simplification  is  made.   (1  +  «~)  ' 

oX 

is  expanded  binomially  to  yield: 

(2.10) 

The  first  two  terms  only  are  retained  since  this  factor  multi- 
plies the  right  hand  side.   The  resulting  equation  is: 


3Z£     1   3  £  +  _vb   a^|   _  ,ODN  3?  8Z£ 
-o   °t    C0 


Sx"7   CT7"  dt*  +   C^  J^t  ~    (23)Tx"  3x^"  (2*11} 
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where  3  =  (y+l)/2  for  gases  and  3=1+  B/2A  for  liquids.  This 
is  in  essence  the  equation  treated  by  Keck  and  Beyer  [4 ] ,  a- 
long  with  the  boundary  condition  £(0,t)  =  -50coswt,  in  their 
perturbation  analysis,  after  terms  of  0((a/k)^)are  eliminated. 
The  Keek-Beyer  solution  is  computationally  simple,  but  it  is 
severely  limited  by  the  labor  required  to  calculate  the  terms 
of  the  perturbation  series.   In  fact,  according  to  the  auth- 
ors, the  perturbation  series  is  probably  only  convergent  for 
3ek(l-e     )/4a  <  1,  which  restriction  corresponds,  for  exam- 
ple, to  an  acoustic  pressure  of  .125  atmospheres  at  1  MHz  in 
water;  one  could  easily  argue  that  this  pressure  hardly  qual- 
ifies as  finite  amplitude, but  then  if  one  is  able  to  measure 
finely  enough,  any  acoustic  signal  can  be  considered  to  be  of 
finite  amplitude.   Actually,  the  series  has  been  demonstrated 
to  be  everywhere  convergent  but  is  still  limited  in  applica- 
bility, by  truncation,  to  values  of  the  parameters  as  defined 
above  [7].   It  is  to  be  noted  also  that  3ek/4a  =  T/4 ;  thus  the 
Keek-Beyer  solution  will  not  allow  T  >  4,  and  since  F  >  4.5  is 
required  for  shock  formation  in  the  viscous  plane  wave  case, 
this  solution  deals  with  non-shocked  waves  only. 

Blackstock  has  demonstrated  that  the  governing  equation 
may  be  reduced  to  a  Burgers   equation  by  applying  the  same 
assumption  utilized  in  the  perturbation  analysis.   This  as- 
sumption, in  fact,  is  the  same  one  which  allows  one  to  arrive 
at  an  analytic  steady  state  solution  of  the  following  problem: 
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It*  "  a^Ft  "  co  a^  '  (2*12) 


namely  the  small  signal  problem  for  plane  wave  propagation  in 
a  viscous  fluid.   This  assumption,  that  a/k  <<  1,  is  valid  for 
most  fluids,  with  the  exception  of  the  most  viscous  ones,  at 
any  reasonable  acoustic  frequency.  In  any  case,  the  effect  of 
making  this  assumption  is  to  allow  one  to  utilize  the  expres- 
sion: 


■^  = L  M.  to    13) 

9x      C0  3t  '  iz.xj, 


in  reducing  the  full  equation  to  a  Burgers  equation.   This 
expression  is,  of  course,  the  differential  equation  which  des- 
cribes the  propagation  of  a  small  amplitude  wave  in  a  lossless 
medium.   The  reduction  procedure  follows.   Starting  with  Equa- 
tion 2.11  (an  approximate  one  which  has  been  obtained  by  as- 
suming a  small  Mach  number),  one  first  substitutes  E,      =  ~—  £. 


o 


on  the  right  hand  side  to  obtain 


'xx  ~  C-7"  ^tt  +  C~^  ^xxt      C   Vxx   '  (2.14) 

o         o  o 


Then  note  the  following  identities  from  Equation  2.13: 


?xt  -  -  F-  5tt  (a>        Co?xtt  -  -  «ttt  (C) 

o 

C„Z      4.  =  "  £  4.*.  (b)         C2£   .  =  -  C  %    ..     =  £...   (d) 

o^xxt  ^xtt  o^xxt       o^xtt    ^ttt 
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K        =  -  7T-   £  4.    (e)      K        =  ^r  ^4.4.    (f)     (2.15a-f) 
^xx      C    xt  xx    C  *      tt 

o  o 


Substituting  (d)  into  the  viscous  term  of  Equation  2.14, 


vb     r  vb       lr  ,  vb 


(£_,,.    +    S„„J    =    o£-£-    (£...     -   C    E   ^J    and 


C    z    Sxt         2C    z     vq,xxt         ^xxt;  2C    *     ^ttt         "o^xtt 


^xx   ~    C^tt   +   55-T    ^ttt    "    Co^xtt)    =    "   C~  Mxx       (2-16) 
o  o  o 


Substituting  (e)  and  (f)  into  the  non-linear  term  of  Equation 
2.16, 


C   HSx    C   v  H^xx    HSx' 
o  o 


C   *C   ^t^xt    C-2"  Vtt* 
o   o  o 


The  equation  after  substituting  and  manipulating  signs  be- 
comes 

Co^tt   +   f    (Co?xtt   -    W    -    Co^xx   =    *Wtt   -    WxtJ 

(2.17) 

Performing  an  operation  which  is  the  integral  equivalent  of 
differentiating  with  respect  to  the  operator  -t-t-   -  C  -r—  ,  one 

at      O  oX 

obtains 

Coh  -   f  6tt  +  C^x  =  I  BCoCt2  (2.18) 
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Dif ferentiating  with  respect  to  ' t'  and  substituting  u  =  £.  , 


vb       „3 


C^Uo.  ~  ^~   U4.4.  +  C  u   =  3C  u  Ux.  (2.19) 

o  t    2   tt    o  x     ou  ut 


Finally,  with  the  variable  transform 


t'  =  t  -  — 
C 
o 


one  obtains 


3.      ~~         1 


Coux'  "  3Cou  uf  =  2  Vbut't'    '  (2'20) 


which  is  Burgers   equation  for  a  boundary  value  problem.   The 
significant  point  here  is  that  one  is  considering  an  equation 
which  is  non-linear,  and  which  includes  viscous  terms.   The 
approximation  a/k  <<  1,  is  applied  two   times  to  reduce  the 
equation  to  Burgers   form,  each  time  to  simplify  viscous  and 
non-linear  terms.   This  is  indeed  an  acceptable  engineering 
approximation;  however,  it  certainly  is  desirable  to  have  the 
ability  to  solve  the  equation  without  resorting  to  this  ap- 
proximation, particularly  for  highly  viscous  fluids  and  per- 
haps certain  dielectric  crystals  with  large  stress  coeffi- 
cients . 

As  noted  by  Blackstock  [6] ,  the  Burgers   equation  solu- 
tion has  one  further  limitation  which  might  be  considered  more 


-22- 


severe.   The  steady  state  solution  to  the  Burgers  equation 
is 


C  =  Y  e  (-l)ni  (i  De"n  0/rcos  ny  (2.21) 

^,  nn2  J 

n=l 


where  L,    is  related  to  u  by, 


u  =  UQ  ~  [log  c]  (2.22) 


Here,  e   is  the  Neumann  factor  (e   =1,  all  other  e   =  2) ;  I 
'   n  o  n        n 

is  the  nth  order  Bessel  of  imaginary  argument;  a  =  3ekx,  a 
non-dimensional  space  variable;  and  y  =  wt',  a  non-dimensional 
temporal  variable.   For  values  of  r  >  50,  the  series  represen- 
tation of  the  solution  is  very  slowly  convergent,  and  there- 
fore difficult  to  use  practically  in  numerical  analysis.   Con- 
sequently, for  T  >  50,  Blackstock  utilized  an  asymptotic  ap- 
proach to  the  solution,  which  in  fact  is  not  valid  near  the 
origin.   To  solve  the  Burgers  equation  for  Bekx  <<  V,    he  re- 
sorts to  a  perturbation  analysis  in  Y      ,  but  recalling  V  =    3ek/a, 
this  is  in  some  sense  similar  to  the  Keek-Beyer  approach  which 
utilized  e  as  the  perturbation  quantity  to  solve  the  same  basic 
equation.   Blackstock 's  solution  for  small  x  contains  terms  of 
first  order  in  the  perturbation  parameter;  for  Y    >>    o ,    this 
should  indeed  be  appropriate,  however  the  coefficients  derived 
contain  an  infinite  series  of  Bessel  functions  which  may  be 
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slowly  convergent  themselves.   Herein  lies  another  motivation 
for  the  parametric  differentiation  solution,  the  desire  to 
produce  a  unified  solution  applying  the  same  assumptions  and 
having  the  same  limitations  throughout. 

In  any  event,  there  is  no  simple  transform  which  will  re- 
duce either  (1)  the  full  equation  valid  for  high  Mach  number, 
or  (2)  cylindrical  and  spherical  finite  amplitude  wave  equa- 
tions to  the  Burgers   form.   Therefore,  an  alternative  ap- 
proach to  the  problem  is  required. 
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III.  PARAMETRIC  DIFFERENTIATION 

Parametric  differentiation  is  a  procedure  which  allows  a 
non-linear  equation  which  involves  a  parameter  to  be  solved  by 
transforming  the  equation  to  a  linear  equation  in  parameter 
space.   This  is  particularly  valuable  for  cases  in  which  nu- 
merical solution  will  be  required  because  it  eliminates  the 
possibility  of  obtaining  multiple  roots  in  solving  non-linear 
difference  equations.   It  is  essential  to  note  at  this  point 
that  some  insight  into  the  nature  of  the  solution  at  the  ex- 
tended range  of  the  parameter  of  interest  would  be  extremely 
helpful.   If  the  solution  is  singular  in  the  parameter,  that 
is,  if  the  nature  of  the  solution  changes  radically  as  one  ap- 
proaches the  limiting  value  of  the  parameter  in  the  known  base 
solution,  then  the  technique  is  not  applicable.   Thus,  para- 
metric differentiation  would  never  be  applicable  to  the  entire 
class  of  singular  perturbation  problems.   The  method  of  appli- 
cation follows. 

Suppose  one  is  given  the  general  problem: 

N[(|)Cx,  t;  i>)]    =   0  (3.1) 


<j>[  (XB,  t;  ij;)]  =  6B  (t;  i>)  (3.1a) 


(J)  1(XB,  t;  x\))]    =  6B  (t;  ty)  (3.1b) 


a<t>[(xB,  t;  4>)]    +  C(j>n(xB,  t;  \p)    =  0B  (t;  ty  )  (3.1c) 


-25- 


where  N  is  a  non-linear  operator;  X  is  a  spacial  vector  coor- 
dinate; t  is  a  time  coordinate;  if;  is  the  parameter  of  inter- 
est; and  the  6  '  s  are  functional  values  of  $> ,    <b    ,    or  a6  +  cd> 
B  nn 

at  XD  as  appropriate  to  the  problem,  n  indicating  the  normal 

13 

derivative.   Initial  values  could  also  be  prescribed  if  appro- 
priate.  If  a  solution  cf>  =   <j>(X,t;  if;  )  exists  at  some  limiting 
value  of  the  parameter  if;  ,  and  it  is  desired  to  obtain  a  solu- 
tion to  the  same  problem  with  identical  boundary  and/or  ini- 
tial values  at  if;   +  A  if;,  then  parametric  differentiation  may  be 
applied.   First,  differentiate  the  governing  equation  and 
boundary  and/or  initial  conditions  with  respect  to  the  para- 
meter to  yield: 


L[g(X,  t;  if;)]  =  0  (3.2) 


g(>tB,  t;  if;)  =  gfi(t;  if;)  (3.2a) 


g  (X  ,  t;  \\>)    =    g   (t;  if;)  (3.2b) 

n   B  B2 


ag(XB,  t;  if;)  +  c9n(xB'  t>    *)  =  9B  (t;  i|»)  (3.2c) 


where 


g(X,  t;  if;)  =  g-  (<}>(X,  t;  if;))  (3.3) 

One  should  realize  that  this  procedure  cannot  be  applied  if 
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the  parameter  appears  in  the  basic  equation  raised  to  a  power 
less  than  one,  and  if  the  limiting  value  of  the  parameter  \\j 
=   0. 

The  operator,  L,  is  a  new  linear  operator  obtained  through 
differentiation.   The  new  equation  with  associated  boundary 
and/or  initial  conditions  may  now  be  solved  in  parameter  space, 
either  analytically  or  numerically,  for  successive  values  of 
the  parameter  \p .      As  the  equation  is  solved  for  each  value  i> , 
a  quadrature  is  performed  which  solves  Equation  3.3  and  re- 
turns a  value  or  functional  form  which,  when  applied  to  the 
solution  for  ty  <  _■,  r    yields  the  desired  solution  at  ip .  .   There 
is  no  numerical  restriction  on  the  range  of  \p;    however,  there 
certainly  may  be  a  physical  limit,  and  some  insight  is  re- 
quired here  to  avoid  blind  application  of  the  method  to  obtain, 
for  instance,  solutions  for  a  range  of  \\)   which  is  beyond  the 
region  of  validity  of  the  basic  equation.   The  quadrature  may 
be  expressed  as  follows: 

fi 

*±  (X,  t;  i() i)  =  <f>i_;L  (X,  t;  *i_1)  +     gcji|)  (3.4) 

♦i-1 

For  more  complicated  equations  which  require  numerical  solu- 
tion, this  step  can  be  quite  complex,  depending  on  the  form  of 
the  g  versus  \p   curve.   In  some  cases,  trapezoidal  integration 
may  be  sufficient;  in  others,  higher  order  integration  techni- 
ques may  be  required.   It  is  clear  in  examining  the  form  of 
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the  equation  above,  that  this  is,  in  essence,  a  recursion  type 
relation  and  consequently  it  demands  a  known  initializing 
function  for  i=l. 

In  this  work,  the  initializing  function  is  the  Fubini 
solution  to  the  problem, 


<1  +  i»26  P  =  co0  <3-5' 


with  the  boundary  condition  £(0,  t)  =  -  H  coswt,  this  being  a 
duplicate  of  Equation  2.9,  without  the  dissipative  term.  Thus, 
the  parameter  of  interest  here  will  be  some  non-dimensional 
quantity  which  contains  v,  and  the  quadrature  will  be  per- 
formed with  respect  to  the  non-dimensional  quantity.   An  ex- 
plicit form  of  the  Fubini  solution   was  given  by  Keck  and 
Beyer  [4]  in  terms  of  particle  velocity, 

00     J    (ni<x) 
u(x,    t)    =    211      I      -J1-~- —  sinn(cot   -   kx)  (3.6) 

n=l 


where  J   is  the  nth  order  Bessel  function  and  k  =  3ek.   This 
n 

solution  is  limited  in  Mach  number  (p^°-  <<  1)  since  a  binomial 

co 

expansion  in  Mach  number  was  utilized  to  arrive  at  this  result. 
An  exact  implicit  solution  is  available  which  could  be  used  for 
analysis  at  greater  Mach  numbers, 

_(I±1) 

u(x,    t)    =   uosin(tot   -  ijp    (l   +  XZ±.  i.)       y    X    )  (3.7) 

o  o 


-28- 


However,  the  numerical  problems  associated  with  an  iterative 
solution  at  many  points  would  be  difficult  to  resolve.   This 
solution  predicts  the  formation  of  the  shock  at  the  point 
where  the  velocity  profile  becomes  negatively  infinite,  that 
is,  where  8u/3x  ->  -°°.   Differentiating  Equation  3.7  with  res- 
pect to  x, 

Y+l  Y+l 

o  o  o  o 


y+l    wx    M         y-1    u   v       Y-l  Y-l    3u    ,  , .    R, 

+  :Fic"(~2~cT)  2C~^]  (3'8) 


or 

||  =  _%»£    tl   +lTi  JL)      Y-l  cos  (      )/(l-lii^ 

O  O  O        O 


-  (3C±i)  -i 

(1  +  V1  i~}    Y  x     cos(    )  ]  (3-9) 


This  quantity  becomes  infinite  when 

2Y 

x=  (i  +  lzi  u)Y-1/1+1Sq^-Co8(   )  (3.10) 

o  o  o 

Noting  that  3  =  j~i    £  =  t^2-/  tt-   =  k,  and  that  the  first  shock 

o   o 

will  occur  at  the  head  of  the  waveform  where  cos (   )  =  1,  and 
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also  that  u  =  0  at  the  shock,  the  discontinuity  distance  X 
=  tt — —  is  obtained.   X  is  the  limit  of  convergence  of  the  ser- 
ies  solution  stated  above.   Thus,  using  the  Fubini  explicit 
solution  as  a  starting  solution,  one  is  limited  to  application 
of  the  technique  of  parametric  differentiation  in  the  region 
0  <   x  <  x  and  for  Mach  numbers  which  are  small  compared  to 
unity.   It  is  clear  from  this  discussion  that  one  is  fundamen- 
tally limited  when  applying  parametric  differentiation  by  the 
quality  of  the  base  solution.   It  is  unfortunate  that  in  this 
case,  the  base  solution  is  limited  with  respect  to  Mach  number; 
on  the  other  hand,  signals  of  sufficiently  large  amplitude  may 
no  longer  be  thought  of  as  acoustic  waves  and  may  propagate  in 
a  different  manner. 


-30- 


IV.   APPLICATION  OF  THE  METHOD 

It  is  one  thing  to  discuss  the  use  of  parametric  differ- 
entiation in  abstract  terms,  and  quite  another  to  execute  its 
application,  there  being  no  real  'theory1  with  respect  to  how 
it  should  be  applied.   In  any  case,  the  first  step  in  the  sol- 
ution technique  is  similar  to  that  for  most  other  numerical 
solution  procedures,  that  is,  non-dimensionalize  the  equation. 
Starting  with  the  reduced  equation 


«xx  -  C^  ?tt  +  ^   W  =  2BVxx  ,  (4.1) 

o         o 


one  introduces  characteristic  lengths  and  times.   This  being  a 

wave  propagation  problem, it  is  reasonable  to  select  x   =  A  and 

t   =  —  ,  where  A  is  the  wavelength;  then  the  non-dimensional 

c   Co 
variables  are: 


x  =  x/x   =  x/A     t  =  t/t   =  tC  /A    and   £  =  £/x   =  Uk 
c  c     o  w  c    ^ 


Then,  Equation  4.1,  in  terms  of  dimensionless  variables,  is 


£**  -    £:£  +  5TT-  E^£  =  23CS"  (4.2) 

^xx    ^tt   C  A  sxxt      sxsxx 


This  equation  is  as  one  would  expect  for  a  'good'  non-dimen- 
sionalization;  the  coefficients  of  the  first  two  terms,  which 
in  fact  should  dominate,  are  of  0(1).   Alternatively,  with  the 
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equation  in  the  form 


(1  +  \    )2^{£..     -  vb£   J  =  C2£  (4.3) 

v     sx     ^tt      bxxt     o^xx 


a  similar  non-dimensionalization  will  yield 


*«    9  ft  ^        \<V\       A  A 

f1  +  ej>  B«tt  -cT:«™t»  "  ^  (4-4) 


which  demonstrates  that  the  relatively  large  numerical  value 
of  the  coefficient  of  the  non-linear  terms  in  4.2  is  simply  a  re- 
sult of  the  expansion  performed,  and  should  not  be  considered  im- 
proper.   The  viscosity  coefficient  appears  as  ^— r-  which  is 
like  an  inverse  Reynold's  number.   However,  because  C   is  a 
propagation  speed  and  not  a  particle  velocity,  this  is  not  a 
classic  acoustic  Reynold's  number.   In  any  case,  this  quantity, 
which  will  be  termed  '^',  is  the  parameter  of  interest,  on 
which  the  quadrature  will  be  performed.   The  solution  for  ^=0 
is  known;  it  is  the  Fubini  solution.   Then  one  has 

oo 

u    (x,    t;    ib    )    =   2U      T   ^n{nKx)  sinn(a3t   -    kx)  (4.5a) 

o  ro  o    L.      ni<x 

n=l 

or,  in  terms  of  particle  displacement 


£  (X,  t;  i[>  )  =  -2^°-  I    Jn(,n>cx)cosn(cot  -  kx)  (4.5b) 

n=l 

which  corresponds  to  the  function  <J)   in  the  theoretical  devel- 
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opment  of  parametric  differentiation,  and  with  5   =  — — . 

The  next  step  in  the  development  is  to  differentiate  E- 
quation  4.2  with  respect  to  ty, yielding 

g   -  g^.4.  +  ^g  4.  +  C  x.  =  23^  g   +  23?  g        (4.6) 

^xx    ^tt    r^xxt    xxt      x^xx      xx^x 


where 


g-||  (4.7) 


and  the  umlaut  notation  has  been  dropped  from  the  non-dimen- 
sional variables.   One  would  immediately  look  at  this  equation 
and  express  horror  at  what  has  been  done;  there  seems  to  be 
added  complexity  rather  than  simplification.   However,  one  can 
now  treat  the  equation  as  an  inhomogeneous,  linear  equation 
with  variable  coefficients,  which  should  make  it  more  tract- 
able for  numerical  solution.   The  only  non-linearity  which  re- 
mains is  found  in  Equation  4.7.   It  should  be  noted  that  the 
non-dimensionalization  in  terms  of  the  size  of  the  coefficients 
has  not  been  adversely  affected.   In  fact,  the  dominance  of  the 
first  two  terms  has  been  enhanced  since  £  ,  £   ~£  and  e  <<  1. 

X     XX 

Alternatively,  differentiating  Equation  4.4  one  obtains 


(2BX1  +ex)2B_1gx(5tt  -  *5xxt>  +  d  +  Cx)26(gtt 


"  *^xxt  "  W  "  9xx  (4-8> 
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which  really  looks  as  if  this  "simplification"  could  have  been 
done  without.  However,  if  one  multiplies  through  by  (1  +  £  ) , 
the  following  equation  is  obtained 


(23)  (1  +  5X)  23(^tt  "  ^xxt)gx  +  (1  +  y  23+1<9tt  -  **xxt 


-£.)=(l+£)g         or 

xxt         ^x  ^xx 


2^xA  +  (1  +  ^x^23+1^tt  -  *W  -  W  =  (1  +  SXK 


(4.9) 


which  is  of  the  same  form  as  Equation  4.6.   This  now  is  the 
equation  which  will  be  utilized  to  solve  the  finite  amplitude 
plane  wave  propagation  problem  without  approximation,  save 
those  which  may  be  necessary  to  obtain  a  base  solution.   In 
this  case,  this  allows  a  significant  extension  because  this 
equation  does  not  require  the  approximation  a/k  <<  1,  and 
therefore,  the  Fubini  solution  may  be  used  as  a  base  solution 
for  the  problem  of  propagation  in  highly  viscous  fluids  since 
it  is  limited  by  Mach  number  only. 

The  quadrature  required  to  extract  the  desired  solution 
is  essentially  dependent  on  the  problem  itself.   Landahl  and 
Ruppert  [13]  employed  Runge-Kutta  integration  to  solve  a  boun- 
dary-layer problem.   In  many  cases,  it  is  not  necessary  to  re- 
sort to  such  exotic  techniques;  in  this  analysis,  g  varies  quite 
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slowly  for  small  values  of  i>;    a  rapid  non-linear  change  in  g 
with  \\)   only  occurs  as  \p   becomes  very  large,  that  is,  as  the 
viscosity  becomes  high.     For  most  fluids,  including  very 
viscous  ones,  a  value  of  ip  sufficient  to  produce  rapid  varia- 
tion of  g  with  \\>   would  not  be  attained  at  any  reasonable  acou- 
stic frequency.   Consequently,  the  integration  technique  em- 
ployed in  this  analysis  was  a  very  simple  straight  line  inte- 
gration, which  required  one  iterative  step.   The  initial  sol- 
ution of  the  equation  for  g  on  the  ith  step  was  computed  based 
on  the  ith  value  of  \\)   and  a  solution  E,  .  (x ,  t)  which  was  ob- 
tained in  the  following  manner 


K±   (x,  t)  =  ^i_1(x,   t)  +  gi_1(x,  t)  •  Uk  -  ^±-1)   (4.10) 


With  this  new  solution,  g. (x,  t) ,  a  modification  of  £ . (X/  t) 
was  performed  as  follows: 


Z±    (x,  t)  =  £±    (x,  t)  +  |(gi(x/  t)  -  gi_1(x,  t)) 


•  (*±  -  ipi-:L)  (4.11) 


The  net  result  being 


C±(x,  t)  =  Ci_1(x,  t)  +  |(gi(x,  t)  +  gi_1(x/  t)) 


(*±  "  ^x)  (4.12) 
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which  is  in  essence  a  trapezoidal  integration.   Clearly,  this 
scheme  would  not  be  useful  for  g  functions  which  vary  rapidly 
and  non-linearly  with  ip ,  that  is,  if  Equation  4.7  exhibits 
highly  non-linear  behavior.   This  scheme  could  be  used  in  an 
iterative  manner,  however,  with  continuing  refinement  of  the  g 
values  at  each  step.   The  entire  process  used  to  arrive  at  the 
extended  solution  for  £(x,t)  is  shown  schematically  in  Figure 
2.   It  may  be  noted  that  there  is  a  requirement  to  obtain  var- 
iable coefficients  at  each  step,  this  of  course  is  a  situation 
which  is  not  enviable  in  a  numerical  solution  in  which  the  co- 
efficients cannot  be  described  analytically,  particularly  when 
these  variable  coefficients  involve  derivatives  of  the  depen- 
dent variable.   In  this  case,  all  the  coefficients  involve  de- 
rivatives. A  simple  solution  to  this  problem  has  not  been 
found;  however,  in  this  analysis,  numerical  differentiation  of 
the  E,    solution  to  obtain  the  coefficients  does  not  seem  to 
have  affected  the  results. 

Another  not  so  simple  problem  in  the  application  of  para- 
metric differentiation  is  the  treatment  of  the  boundary  condi- 
tions.  One  must  be  extremely  careful  in  not  trying  to  read 
too  much  into  the  problem.   In  this  analysis,  two  approaches 
were  considered  for  determination  of  boundary  values  for  appli- 
cation in  the  finite  difference  scheme,   the  result  being  that 
the  most  natural  and  straightforward  approach  was  eventually 
chosen  as  the  proper  one.  In  analyzing  this  problem,  one  can 
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make  some  subjective  evaluations  of  the  expected  behavior  of 
g  based  on  what  the  gAi|>  product  represents.   Dissipation  will 
tend  to  slow  the  non-linear  generation  of  harmonics  as  the 
wave  propagates,  and  if  dissipation  is  strong  enough,  one 
would  expect  the  waveform  to  remain  sinusoidal,  if  such  is  the 
nature  of  source  excitation,  no  matter  how  strong  the  non-lin- 
earity is  in  absolute  terms,  that  is,  for  T  <  1  the  waveform 
should  not  steepen  appreciably.   In  any  case,  near  the  bound- 
ary the  effect  of  gAiJj  should  be  to  decrease  the  amplitude  of 
the  signal.   Aip  is  always  positive,  therefore  to  obtain  such 
a  reduction  in  amplitude,  one  would  expect  the  sign  of  g  to  be 
opposite  that  of  the  base  solution.   Attempting  to  produce  an 
analytical  model  which  fulfills  one's  expectations  near  the 
boundary  may  be  self  defeating.   In  this  analysis  it  was  ini- 
tially considered  desirable  to  introduce  a  small  perturbation 
of  the  boundary  condition  in  terms  of  ty,    and  consequently  the 
boundary  condition  was  rewritten  as, 


HO,  t)  =  -  (~°   cos(a)t)  (4.13) 


which  yielded  the  intuitively  desirable  result  in  terms  of 
sign, 


^(UO,  t)  )  =  g(0,  t)  =  {1l^yz    cos(o)t)  (4.14) 


This  model  also  seemed  appropriate  since  it  produced  a  bound- 
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ary  value  of  g  which  decreased  slightly  as  \p    increased,  corres- 
ponding to  an  anticipated  increase  in  the  viscous  effects.   If 
one  recalls  the  nature  of  the  quadrature,  the  only  effect  of 
this  model  would  be  to  force  the  solution  near  the  boundary  to 
take  on  the  desired  form.   This  'forcing1  of  the  solution  by 
injecting  one's  subjective  expectations  into  the  modelling  of 
the  boundary  conditions  was  not  particularly  bad,  but  neither 
did  it  accomplish  anything  with  respect  to  improving  the  solu- 
tion.   If  one  instead  blindly  applies  the  tt-j-  operator  to  the 
boundary  condition,  the  result  is 


g(0,  t)  =  ^(-Hocos(ojt)  )  =  0  (4.15) 


which  does  absolutely  nothing  to  convince  one  that  the  proper 
sort  of  quadrature  will  occur  near  the  boundary,  but  instead  re- 
lies on  the  differential  equation  to  properly  predict  the  g 
values  everywhere.   In  the  final  analysis,  this  was  the  cor- 
rect procedure;  there  were  in  fact  very  small  differences  in 
the  solution  values  obtained  with  the  forced  boundary  condi- 
tions, and  those  obtained  with  the  'natural'  condition. 

One  may  wonder  why  initial  conditions  are  not  important 
in  this  problem  since  it  is  essentially  hyperbolic  in  char- 
acter, and  such  problems  are  generally  thought  of  as  initial 
value  problems.   The  answer  is  that  when  one  is  interested  in 
the  steady  state  solution  only,  for  a  problem  with  a  periodi- 
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cally  varying  boundary  condition,  and  in  this  case  that  is  the 
solution  desired  for  £ (x ,  t) ,  one  may  set  the  condition  that 
the  initial  values  occurred  at  some  sufficiently  large  time  in 
the  past  that  all  start-up  transients  have  disappeared,  and 
therefore  only  the  boundary  conditions  determine  the  behavior 
of  the  solution  [16].  It  is  only  natural  to  expect  this  physical 
interpretation  to  be  carried  over  into  the  parameter  space  ap- 
plication, and  consequently,  initial  values  have  been  of  no 
concern. 

It  is  valuable  to  note  that  the  final  result  for  the  cho- 
sen maximum  value  of  ty   need  not  be  the  only  output.   For  a 
given  $£  product,  a  solution  for  any  r  consistent  with  the  ac- 
curacy obtainable  (as  r  -*-  °°,  the  accuracy  required  to  obtain 
meaningful  results  is  increased  for  a  given  Be  because  the 
differences  between  the  inviscid  and  viscous  solutions  go  to 
zero)  ,  may  be  obtained  by  selecting  the  appropriate  ty   -  ^-jr   and 
simply  calling  for  output  at  the  end  of  the  second  iterative 
step.   For  a  given  fluid  this  would  allow  one  to  analyze  the 
losses  at  different  frequencies  for  a  given  Be  and,  in  this 
case,  out  to  the  discontinuity  distance.   In  fact,  tabulations 
of  g  versus  ip  at  specific  points  could  be  stored  so  that  the 
solution  could  be  extracted  in  the  future  simply  by  integra- 
ting and  applying  the  correction  at  the  appropriate  point  in 
the  base  solution. 

It  is  useful  to  compare  the  base  solution  waveform  with 
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the  anticipated  extended  solution  waveform,  the  comparison  in 
this  case  being  very  simple  because  the  extended  solution  or 
a  very  good  approximation  thereto  is   well  known.   As  it  de- 
velops, the  incremental  differences  between  the  base  solution 
and  the  extended  solution  are  expected  to  be  in  phase  with  one 
another,  if  they  are  written  in  the  form 


A£(x,  t)  =  Cf(x,  t)  -  CQ(x,  t)  (4.16) 


where  A£(x,  t)  is  the  difference  between  the  calculated  solu- 
tion and  the  base  solution  when  integrating  Equation  4.7  from 
ip   to  ijjf.   The  observance  of  A£(x,  t)  varying  in  the  expected 
manner  is  a  key  clue  to  the  success  of  the  application  of  the 
method.   This  behavior  was  in  fact  observed  in  this  analysis, 
and  will  be  discussed  in  a  following  section.   Once  one  has 
developed  an  equation  in  parameter  space  and  made  a  subjective 
evaluation  of  what  he  expects  to  be  the  functional  form  of  g 
and  of  A£(x,  t) ,  he  is  ready  to  proceed  to  the  numerical  tech- 
nique to  be  used  for  solution.   Ideally,  it  would  be  hoped 
that  the  transform  to  parameter  space  would  lead  to  an  anal- 
ytic solution  of  Equations  4.6  or  4.9  and  Equation  4.7,  or  at 
the  very  least,  Equation  4.6  or  4.9.   Such  is  not  the  case 
here,  and  numerical  solution  of  both  was  required.   Thus,  to 
analyze  the  results,  and  prove  the  worth  of  the  method,  a  sub- 
jective feel  for  the  results  was  required. 
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V.    FINITE  DIFFERENCE  EQUATIONS 

After  the  equation  of  motion  has  been  derived,  non-dimen- 
sionalized  and  its  analog  in  parameter  space  has  been  obtained 
through  differentiation,  it  may  still  be  a  significant  task  to 
obtain  the  solution  of  the  linear  g  equation.   In  this  case, 
Equation  4.6  or  4.9  describes   the  function  g  with  boundary 
condition  g(0,  t)  =  0.   In  analyzing  the  equations,  it  is  noted 
that  they  are  essentially  wave-like  in  nature,  at  least  close 
to  the  boundary,  since  the  g    and  g    terms  dominate.   Adopt- 
ing  the  approximation  E,     ~  £.  (in  non-dimensional  form)  and 

.A.         L- 

differentiating  with  respect  to  x,  one  notes  that  E,       »  -u 

XX        X 

and  since  u   ->  -°°  at  the  shock  in  the  inviscid  case,  then  E 

x  xx 

+  »  at  the  shock.   This  wave-like  behavior  may  not  then  be  so 
readily  observable  away  from  the  origin  since  the  coefficients 
involving  E,        may  become  very  large  in  the  regions  in  which 
the  shock  is  forming,  or  has  formed  in  an  extended  solution, 
in  which  case  the  g   term  becomes  important.   The  g     term 
remains  small  since  iJj  is  generally  much  less  than  one.   The  g 

X 

term,  of  course,  is  the  one  which  results  from  the  inclusion 
of  non-linear  terms,  and  it  is  entirely  consistent  that  it  be- 
comes important  in  the  regions  of  the  waveform  which  exhibit 
the  result  of  non-linear  propagation,  namely  in  the  region  at 
the  head  of  each  compressive  portion  of  the  waveform  and  at 
the  tail  of  the  raref active  portion.   Once  a  shock  has  formed, 

it  is  conceivable  that  the  g   term  will  alter  the  character  of 

^x 
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the  equation  to  the  extent  that  it  is  no  longer  wavelike;  this 
is  not  of  immediate  concern  since  this  analysis  is  restricted 
to  the  first  region  of  propagation.   Therefore,  in  attempting 
to  solve  this  equation,  it  was  treated  as  if  it  were  wavelike 
throughout.   As  previously  noted,  the  literature  is  rich  with 
various  methods  of  solution  of  the  fundamental  equation  of  mo- 
tion, techniques  which,  if  applied  judiciously  and  with  essen- 
tially the  same  approximations  as  those  utilized  to  solve  the 
equation  of  motion,  could  possibly  be  applied  to  Equation  4.6 
or  4.9  to  yield  analytic  solutions.   However,  because  of  the 
form  of  the  variable  coefficients,  this  analytic  solution 
would  at  best  be  extremely  difficult  and  most  likely  would  be 
impossible;  so  in  this  analysis  solution,  by  finite  differences, 
of  the  parameter  space  equation  was  accomplished,  and  numeri- 
cal integration  of  Equation  4.7  was  performed.   Since  the  e- 
quation  is  essentially  hyperbolic  in  nature,  a  finite  differ- 
ence scheme  appropriate  to  a  hyperbolic  system  was  the  obvious 
choice. 

There  are  two  broad  categories  of  finite  difference 
schemes  for  partial  differential  equations,  which  are  termed 
explicit  and  implicit.   Adopting  a  standard  notation  for  a 
nine  point  difference  mesh,  Figure  3,  an  explicit  scheme,  is 
one  which  has  only  one  unknown  in  the  i+1  column  or  the  j+1 
row,  depending  on  whether  one  is  "marching"  in  the  i  or  j  dir- 
ection.  An  implicit  scheme  is  one  for  which  there  are  two  or 
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more  unknowns  in  the  i+1  column  or  j+1  row.   To  solve  the  ex- 
plicit scheme  for  a  single  value  in  the  future  (marching  dir- 
ection) one  need  know  only  the  values  of  the  dependent  varia- 
bles at  all  present  and  past  positions  on  the  mesh.   The  terms 
future,  present  and  past  apply  in  a  spatial  as  well  as  tempor- 
al sense.   Thus,  for  the  explicit  scheme,  if  one  has  initial 
conditions  valid  for  all  space,  or  if  one  has  initial  condi- 
tions for  a  half  space  and  a  boundary  condition,  a  complete 
solution  can  theoretically  be  generated  numerically  for  a  pro- 
perly-posed problem.   The  implicit  scheme  requires  solution  at 
two  or  more  points  based  on  present  and  past  points.   Thus, 
one  generates  two  unknowns  in  the  first  application  of  the 
mesh,  and  an  additional  unknown  with  each  subsequent  applica- 
tion, yielding  N+l  unknowns  with  N  equations.   Thus,  an  addi- 
tional condition  is  required  in  order  to  apply  an  implicit 
scheme;  generally  it  is  a  boundary  condition  at  some  point  in 
space.   In  any  case,  this  requirement  essentially  restricts 
the  use  of  implicit  schemes  in  infinite  domain  problems  such 
as  this  one. 

Having  been  limited  to  the  use  of  explicit  solution 
schemes,  there  exists  one  additional  problem  of  significant 
import;  that  is,  stability  of  the  scheme.   Unfortunately,  it 
is  difficult  to  consider  stability  until  one  has  actually  de- 
veloped the  finite  difference  equation,  stability  of  the  scheme 
being  a  function  of  the  way  the  scheme  itself  is  developed. 
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Proceeding  to  finite  difference  approximations  of  the  Equation 
4.6  or  4.9,  one  finds  a  literature  which  is  rich  in  references 
to  the  wave  equation  for  infinitesimal  amplitude  waves,  in  a 
lossless  medium,  but  which  has  very  little  in  the  way  of  ref- 
erence to  non-linear  equations  or  equations  for  a  viscous 
medium.  Recalling  that  the  parameter  space  equation  is  essen- 
tially wavelike  as  long  as  the  coefficient  of  the  g   term  is 
relatively  small,  one  can  use  the  finite  difference  approxima- 
tions to  the  wave  equation  as  a  starting  point.   Since  an  es- 
sential aim  of  this  study  was  to  demonstrate  the  utility  of 
the  method  of  parametric  differentiation  in  non-linear  wave 
problems,  a  second  order  finite  difference  scheme  was  selected; 
it  is  clear  that  greater  accuracy  could  have  been  obtained  with 
a  higher  order  scheme.   The  following  finite  difference  approx- 
imations were  used  for  each  of  the  derivatives  of  g  in  either 
Equation  4.6  or  4.9: 


gxx  =  (gi-l,j  "  2gi,j  +  9i+l,j>/(Ax)   +  0((Ax)  >     (5-1} 


*tt  =  (gi,j-l  "  2gi,j  +  9i,j+l)/(At)2  +  0((At)2)     (5.2) 


gx  =  (gi+i,j  gi-],y 


2Ax  +  0( (Ax) 2)  (5.3) 


g     =  (g.     .  -  2g.   .  +  g.,,   .  -  g.  ,   .  n  +  2g .  .  . 
^xxt     ^1-1, j      i/D     l+l# j     l-l # j~l      i/D~- 


-  gi+i/j_i)/Ax2At  +  0((Ax)2(At))  (5.4) 
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where  Ax  and  At  denote  the  step  sizes  in  space  and  time  res- 
pectively. It  must  be  stressed  that  these  approximations,  par- 
ticularly Equations  5.3  and  5.4,  are  not  the  only  ones  which 
could  be  used.   For  instance,  a  forward  or  backward  difference 
could  have  been  used  in  place  of  the  central  difference  of  E- 
quation  5.3.   This  would,  however,  have  changed  the  order  of 
the  approximations.   The  form  of  Equation  5.4  is  that  of  a 
second  central  difference  in  space  and  a  backward  difference 
in  time,  the  backward  difference  being  necessary  in  order  to 
preserve  the  explicit  nature  of  the  scheme.   A  central  differ- 
ence or  forward  difference  in  time  for  Equation  5.4  would  have 
not  only  destroyed  the  explicit  nature  of  the  scheme  but  would 
also  have  yielded  an  unknown  to  be  determined  at  the  point 
g. ,,  .,,  whose  coefficient,  ^,  would  have  been  much  smaller 
than  those  of  the  other  terms,  and  could  have  led  to  an  in- 
stability, even  if  one  were  using  an  implicit  scheme.   This 
problem  could  be  resolved  in  an  implicit  scheme  by  substitu- 
ting for  Equation  5.1  a  weighted  time  average  of  second  order 
approximations  of  g 

XX 

Representing  the  variable  coefficients  of  g   ,  g   ,  g  , 

XX      T-.  X-      X 

and  g  ,  by  A,,  A.,  A.,  and  A.  respectively,  Equation  4.6  or 
4.9  becomes,  upon  substituting  the  finite  difference  approxi- 
mations 5.1-4: 

-A4-gi_1^j_1/(Ax) 2At  +  (A1/(Ax)2  +  A3/2Ax  +  A4/(Ax)2At) 
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*gi-l,j  +  ("A2/(At)2  +  2A4/(Ax)2At)  .gi^j_1  +  (»2A]L/(Ax) 


+  2A2/(At)2  -  2A4/(Ax)2At).gi^  -  A2  .gi  ^    /  ( At)  2 


-  A 


4*91+1,  j_1/(Ax)2At  +  (A-j/CAx)2  -  A3/2Ax  +  A4/(Ax)2At) 


.g.+1/j  =  F  (5.5) 


where  F  is  the  inhomogeneous  term.   At  this  point,  one  may 

look  at  Equation  5.5  and  note  that  it  has  one  unknown  in  the 

j+1  row  if  one  is  marching  in  space,  and  two  in  the  i+1  column 

if  marching  in  time;  so  that  it  is  well  suited  to  treatment  of 

an  initial  value  problem  explicitly,  but  is  implicit  if  used 

for  a  boundary  value  problem.   In  order  to  apply  the  scheme  in 

an  explicit  manner  to  the  boundary  value  problem,  one  must  make 

an  additional  approximation.   To  solve  for  g(i+l,j),  one  needs 

information  in  the  i-1  and  i  columns,  plus  information  at  the 

point  g. ,,  .  ,.   Because  the  third  order  derivative  is  present, 
e  ^i+l, j-1  ^ 

there  is  no  way  to  avoid  this  problem.   However,  considering 

the  fact  that  the  coefficient  of  the  third  order  derivative  is 

small  compared  to  the  coefficients  of  the  other  terms,  it  would 

not  be  unreasonable  to  presume  that  if  one  had  an  approximate 

value  of  g-.-i  •_->  /  then  a  solution  of  reasonable  accuracy  could 

be  obtained  for  g. ,,  .,  which  would  be  based  on  five  pieces  of 

yi+l,D 

accurate  information  and  a  single  piece  of  approximate  infor- 
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mation.   Considering  Figure  4,  which  depicts  the  finite-dif- 
ference grid  with  the  coefficients  of  the  various  points  in 
the  solution  cube  annotated  thereon,  one  notes  that  the  ap- 
proximation required  at  g .  ,  .  ,  need  be  made  only  once  in 
each  spatial  step.   After  solving  for  g.  ,  .at  the  'base'  of 
a  single  spatial  column,  q. ..  .  ,  is  now  known  for  all  subse- 
quent  applications  of  the  cube  as  one  marches  out  in  time.  The 
approximation  used  to  arrive  at  this  "corner  point"  is  a  Tay- 
lor's series  expansion  in  which  the  derivatives  at  the  point 
about  which  the  expansion  is  taken  are  approximated  by  back- 
ward differences  accurate  to  O(Ax).   Again  greater  accuracy  is 
obtainable,  but  was  not  required  because  of  the  relative  small- 
ness  of  the  coefficient.   The  resulting  expression  for  the 
"corner  point"  is, 


g. ,,  .  ,  =  2.5«g.  .  ,  -  2»g.  ,  .  ,  +  .5«g.  ~  .  , 
yi+l, 3-1        ^i/3~l      ^1-1,3-1       ^1-2,3-1 


(5.6) 


It  must  be  noted  that  this  approximation  is  a  potentially  weak 
link  in  the  method,  due  to  the  fact  that  near  the  head  of  the 
waveform  the  x  derivatives  become  large,  and  a  truncated  Tay- 
lor's series  may  not  be  of  sufficient  accuracy.   In  any  case, 
it  did  not  appear  to  significantly  affect  the  nature  of  the 
results  in  this  analysis. 

In  the  discussion  of  explicit  finite  difference  schemes 
one  must  inevitably  be  concerned  with  the  stability  of  the 
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scheme  as  it  relates  to  step  size.   According  to  the  Courant- 
Friedrichs-Lewy  criterion  for  stability  of  finite  difference 
schemes  [17],  the  ratio  of  the  step  sizes  must  be  such  that 
one  does  not  attempt  to  find  the  solution  at  a  point  which  is 
outside  the  region  of  influence  of  the  known  data,  the  region 
of  influence  being  determined  by  the  characteristics  emanating 
from  the  endpoints  of  the  known  data.   Unfortunately,  this 
problem  is  difficult  to  address  here  because  of  the  nature  of 
the  governing  differential  equation.   Typically,  third-order 
partial  differential  equations  are  not  treated  as  such,  but 
are  instead  reduced  to  some  more  tractable  form;  in  this  case 
the  approach  using  the  Burgers   equation  was  such  a  reduction. 
Consequently,  little  appears  in  the  literature  with  reference 
to  stability  for  third-order  equations.   An  unsuccessful  in- 
vestigation was  conducted  to  attempt  to  arrive  at  the  exact 
characteristics  by  reduction  to  canonical  form.   However,  this 
equation  is  wave-like,  and  one  could  reasonably  expect  the 
characteristics  to  be  similar  to  those  of  the  wave  equation, 
and  that  the  stability  criterion  would  be  similar  in  form. 
Working  with  this  foreknowledge  and  applying  a  modified  Von 
Neumann  [17]  stability  test,  an  exact  stability  criterion  will 
be  established.   The  VonNeumann  test  essentially  consists  of 
examining  possible  exponential  solutions  to  a  finite  differ- 
ence scheme  to  determine  when  they  may  grow  without  bound 
given  finite  initial  or  boundary  conditions.   With  this  know- 
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ledge  one  may  adjust  the  stepsize  so  that  the  scheme  will  be 
stable.   Unfortunately ,  the  VonNeumann  test  is  applicable  to 
finite  difference  equations  with  constant  coefficients;  in 
this  case  the  coefficients  are  variable,  and  the  only  recourse 
is  to  treat  the  coefficients  as  if  they  were  constant,  derive 
the  stability  criterion,  and  then  discuss  the  implications  of 
their  variability  with  respect  to  how  it  affects  the  stability 
criterion. 

The  exponential  solutions  to  the  scheme  which  may  grow 
without  bound  will  have  the  form, 


im6  mX  /t-    ns 

g    =  e    e  (5.7) 

^m,n 


where  the  subscript  notation  m,n  corresponds  to  the  notation 
i,j  used  previously.  A  is  any  real  number  and  0  is  an  arbi- 
trary complex  number.   Thus,  boundary  data  may  be  written, 


inA  /c  ON 

gQ/n  =  e  (5.8) 


and  the  right  hand  side  may  be  considered  as  a  typical  term  in 
a  Fourier  expansion  of  the  boundary  data.   Substituting  Equa- 
tion 5.7  into  5.5,  ignoring  the  inhomogeneous  term  since  it 
will  not  affect  stability,  and  dividing  through  by  e    e    , 
one  obtains 


■A4e  l9e  lX/(Ax)2At+  (A.,/ (Ax)2  +  A3/2Ax  +  A4/(Ax)2At)e  1' 
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+  (-A0/(At)2  +  2A,/(Ax) 2At)e"lX  +  (^A^Ax)2  +    2An/(At)2 


-  2A„/(Ax) 2At)  -  A0elX/(At)2  -  A,el9e  lA/(Ax)2At 


+  (Aj^/CAx)2  -  A3/2Ax  +  A4/(Ax)  2At)e10  =  0  (5.9) 


Rearranging  terms  such  that  they  are  grouped  by  coef f icients , 
the  following  results, 


eiG+e~ie-2        iA+e~iA  2        i6  e~iG 

V6  (Ax)2"  }  "  A2(6  (AGt)2  "  }  "  A3(!L-2a! } 

i6   -i6_9   -,   ~iA 
+  V   (Ax)^~  >  'T'  =  °  <5-10> 


Recognizing  the  various  terms  as  trigonometric  forms,  one  ob- 
tains 

-  (4A][sin2|)/(Ax)  2  +  (4A2sin2-^)  /  ( At)  2  -   iA3sin0/Ax 

-  A. (4sin2&  (l-e"lA)/(Ax) 2At  =  0  (5.11) 


Now  solving  for  the  imaginary  portion, 


-A3sin0/Ax  -  4A4sin2-|  sinX/(Ax)2At  =  0  (5.12) 


which  identity  indicates  that  the  product, 
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9 
AxAt  =  -4A.sin  y  sinA/A_sin9  (5.13) 

should  hold,  but  this  is  clearly  absurd  since  it  indicates 
that  as  Ax  ->  O,  At  may  become  arbitrarily  large  as  well  as 
negative.   Reason  and  foreknowledge  of  the  stability  criterion 
for  the  wave  equation  dictate  that  this  result  has  no  meaning. 
Returning  to  Equation  5.11  and  equating  the  real  parts,  one 
obtains 

(-A1sin2|)/(Ax) 2  +  (A2sin2|)/At2  -  (A4sin2|)  (1-cosA) / ( Ax)2 At 

=  0  (5.14) 


or 


2„-   2^ 


A~ (Ax)  sin 


2D  2  2 


Sin  2    A, (At) z+A.At(l-cosA)  (5.15) 


Now  one  argues  that  to  prevent  the  assumed  exponential  solu- 
tion from  increasing  without  bound,  6  must  have  no  negative 
imaginary  part.   However,  because  of  the  quadratic  form,  if  9 
has  imaginary  roots,  they  will  appear  in  conjugate  pairs; 
therefore  9  must  be  real.   Since  A  was  assumed  arbitrary,  this 
requires  the  following  restriction, 

A  (Ax) 2 

0  <  r= ,  .  ,  .  ,  ,  .   .  ,  ,, r-T-  <  1  (5.16) 

—  A, (At) 2+A.At (1-cosA)  — 
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Now,  noting  that  the  maximum  value  of  the  expression  above  is 
obtained  when  1-cosX  =  0,  it  is  clear  that  a  sufficient  ex- 
pression for  the  above  is, 

A, (Ax) 2 

^A^At)-!1  <5-17> 

which  is  of  precisely  the  same  form  as  the  stability  criterion 

for  the  wave  equation,  in  which  case  A,  =  A„  =  1.   In  this 

2  R+l 
case  A,  and  A„  are  variable  with  A,  =  1  +  £   and  A~  =  (1 +  £  ) 
12  1x2        x 

for  Equation  4.9.   Since  £   is  approximately  equal  to  the  lo- 

X 

cal  acoustic  Mach  number,  which  is  taken  to  have  a  value  of 
0(10   )  in  this  analysis,  the  criterion  is  only  slightly  dif- 
ferent from  that  of  the  wave  equation.   However,  it  is  impor- 
tant to  realize  that  the  ratio  of  step  sizes  cannot  be  unity 
as  one  may  have  anticipated.   It  is  also  comforting  to  note 
that  A,.,  the  coefficient  of  the  g   term,  did  not  enter  into 

-J  X 

the  criterion  in  any  way,  lower  order  terms  generally  having 
little  effect  on  stability.   This  statement,  however,  deserves 
further  comment.   Recalling  the  discussion  of  the  size  of  the 
respective  terms  in  the  parameter  space  equation,  the  g   term 
has  a  coefficient,  £.   ,  which  will  become  very  large  in  the 

XX 

region  of  a  shock.  It  is  inconceivable  that  this  term  will 
not  affect  the  stability  of  the  scheme  in  such  a  case.  The 
point  is  that  when  £    becomes  very  large  the  entire  nature  of 

XX 

the  equation  does  change  with  the  g   term  being  the  critical 

X 
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term.   In  this  instance  the  stability  analysis  is  no  longer 
valid,  and  one  would  be  forced  to  rely  on  some  other  stability 
test.   This  problem  is  not  great  here  because  the  analysis  is 
restricted  to  the  first  zone  of  propagation.   Also,  the  dis- 
cretization of  the  displacement  data  limits  the  maximum  value 
that  £    can  attain  in  the  numerical  scheme.   This  is,  as  yet, 
an  unresolved  numerical  problem;  to  allow  £,  t      to  attain  its 
true  value  in  a  nearly  inviscid  fluid,  an  exceedingly  small 
step  size  in  space  would  be  required.   However,  if  one  consid- 
ers the  actual  size   £,    may  attain  in  a  typical  propagating 
wave,  this  problem  may  be  put  in  a  proper  light.   Blackstock 
[5]  has  defined  a  shock  thickness  for  the  viscous  problem, 


T  =  2Acosh~1/i7A  (5.18) 


where  T  is  the  thickness,  defined  as  the  spatial  distance  from 
trough  to  peak  of  the  waveform,  and  A  =  2(1  +  o)/i\T.      As  an 

example,  consider  the  case  when  Y   =  100  and  a  =   1;  then  T  = 

8        — ] 
t  AA   cosh   5tt  -  .08.   In  this  case,  £    changes  from  a  maximum 
IOOtt  ^xx     ^ 

negative  value  to  a  maximum  positive  value  over  a  change  in  x 

of  .08.   Since  £    ^  -u    .    the  chanqe  in  u  over  this  distance 

xx     x 

is  approximately  twice  the  Mach  number,  which  as  noted  before 

-2 

is  taken  to  be  of  0(10   ).   Then  an  approximate  value  of  £. 

may  be  taken  as  Au/Ax  ts  .25;  therefore,  while  the  g   term  does 
become  important,  it  does  not  dominate.   For  larger  values  of 
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r,  it  will  become  increasingly  important,  will  in  fact  domi- 
nate the  equation,  and  will  require  smaller  step  size  near  a 

shock  to  correctly  estimate  £ 

J  xx 

Another  consideration  in  the  finite  difference  scheme  is 
the  handling  of  boundary  values.   Treated  as  an  initial-value 
problem,  there  would  be  no  problem  in  starting  the  solution, 
since  £  and  £   would  be  everywhere  prescribed  for  t  =  0.   In 
this  case,  only  one  boundary  condition  is  available  for  x  =  0, 
the  other  condition  being  that  £  -*■   0  as  x  •*■   °°,  which  is  of  no 
import  here  except  that  it  required  the  ruling  out  of  solutions 
which  grow  in  space.   Because  two  "columns"  of  information  are 
needed  at  the  boundary  to  start  the  solution  for  g,  another 
condition  is  required.   Recalling  the  discussion  in  the  previ- 
ous chapter  concerning  the  modelling  of  boundary  conditions, 
this  problem  of  requiring  an  extra  condition  could  lead  one  to 
expend  great  effort  in  attempting  to  develop  a  consistent  and 
meaningful  extra  condition,  when  in  fact  if  one  simply  applies 
the  method  of  parametric  differentiation  as  it  is  defined,  the 
second  condition  will  become  apparent.   The  point  is  that  the 
boundary  is  a  real  thing,  not  some  abstraction  in  parameter 
space,  and  consequently  the  source  excitation  has  no  depen- 
dence on  ^  whatsoever;  therefore  one  can  in  fact,  without  wor- 
rying about  what  £   is  at  the  boundary  (which,  incidentally, 
would  overspecify  the  problem) ,  write  another  boundary  condi- 
tion, 
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g  (0,  t)  =  0  (5.19) 


This  completes  the  information  needed  to  march  out  the  solu- 
tion of  the  problem.   It  is  interesting  to  note  that  because  g 
and  g   are  zero  at  the  boundary,  the  only  non-zero  term  in  the 
finite  difference  cube  on  its  first  application  is  the  inhomo- 
geneous  term.   Without  this  term,  there  would  be  no  solution 
because  all  terms  in  the  finite  difference  cube  would  be  zero. 
Again,  this  is  consistent  with  intuitive  thoughts  about  the 
form  of  the  solution  for  g;  it  should  in  fact  be  very  small 
near  the  boundary  so  that  when  the  quadrature  is  performed  to 
calculate  the  particle  displacement,  the  net  change  in  the 
value  of  E,    from  the  inviscid  to  the  viscous  solution  is  small. 
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VI.   RESULTS 

The  practical  realization  of  the  preceding  development  is 
a  computer  program  which  solves  the  non-linear  plane  wave  pro- 
pagation problem  in  a  viscous  fluid  and  provides  output  which 
allows  one  to  readily  compare  the  results  obtained  with  this 
method ,  and  those  obtained  by  Keck  and  Beyer  [4].   Their  solu- 
tion was  chosen  as  the  point  for  comparison  because  it  was 
simple  and  easy  to  calculate,  and  because  this  analysis  had  to 
be  conducted  on  a  shoestring  budget.   Unfortunately,  this  did 
not  allow  good  comparisons  for  the  higher  ranges  of  T   because 
the  Keck/Beyer  solution  is  not  sufficiently  accurate  [6], 
lacking  harmonic  content  beyond  the  sixth  harmonic.   One  can 
argue  though  that  the  results  for  large  V    appear  to  be  consis- 
tent with  one's  expectations,  and  that  the  method  appears  to 
generate  a  plausible  progression  of  results  as  the  value  of  T 
is  decreased  (corresponding  to  an  increase  in  ty) .   An  annotated 
version  of  this  computer  program  appears  in  Appendix  A.   It  is 
worthwhile  to  note  at  this  point  that  because  of  the  hyperbolic 
nature  of  the  problem,  the  domain  of  dependence  for  a  given 
point  becomes  larger  and  larger  the  further  away  from  the 
boundary  that  the  point  is  located.   The  result  is  that  if  one 
desires  a  solution  value  N  discrete  steps  away  from  the  boun- 
dary, N2  +  2N  +  1  points  are  required  in  the  solution  grid. 
This  is  a  fundamental  limitation  of  this  method  as  it  is  pres- 
ently formulated,  and  may  be  difficult  to  improve  upon. 
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In  this  analysis  the  Mach  number  was  taken  to  be  rather 
large,  e  =  O(.Ol),  in  order  to  minimize  N  while  still  being 
able  to  generate  useful  results.   B/A  was  taken  as  6,  which  is 
the  approximate  value  of  that  quantity  for  water;  thus  3=4. 
For  a  given  3e  product,  the  discontinuity  distance,  X,  will 
occur  at  a  point  which  is  1/2tt3£  wavelengths  from  the  boundary, 
regardless  of  frequency.   Thus,  for  3e  =  .04,  X  »  4X.  In  water, 
a  Mach  number  of  .01  corresponds  to  an  acoustic  pressure  given 
by  the  expression  p  =  p  C2e  or  p  —  225  atmospheres,  or  in  deci- 
bel units,  SPL  =  267  dB  re  1  uPa,  which  is  quite  a  strong  sig- 
nal.  The  values  of  ij;  have  been  varied  from  zero  to  .01  which 
would  reflect  a  change  from  inviscid  to  highly  viscous  propa- 
gation.  Recalling  the  form  ty   =  vb/C  A,  it  is  clear  that  by 
fixing  ty,    the  frequency  for  which  the  solution  is  valid  is  al- 
so fixed  for  a  given  fluid.   To  put  these  results  in  a  proper 
perspective,  a  value  of  \p   =  .  01  in  water  would  correspond  to  a 
frequency  of  approximately  five  gigahertz.   However,  the  re- 
sults for  smaller  values  of  \p   represent  more  reasonable  fre- 
quencies.  Another  interpretation  of  the  results  for  \p   =  .01 
is  that  a  solution  has  been  obtained  for  the  fluid  with  an  ex- 
tremely high  viscosity  coefficient  such  as  glycerin,  a  solu- 
tion which  has  not  previously  been  obtained  analytically,  be- 
cause of  the  restriction  a/K   <<  1,  which  is  necessary  to  re- 
duce the  equations  to  analytically  tractable  forms.   One  might 
conclude  from  these  results  that  there  is  nothing  significantly 
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different  in  the  form  of  the  solution  for  such  a  liquid.   In 
any  case,  the  results  which  follow  will  indicate  the  potential 
of  this  method  for  application  over  a  great  range  of  values  of 
the  non-linear  and  viscous  parameters. 

Figures  5a,  b,  c,  and  d  are  plots  of  the  particle  dis- 
placement calculated  by  parametric  differentiation,  the  par- 
ticle velocity  calculated  from  this  solution,  the  difference 
between  the  inviscid  (Fubini)  and  the  viscous  (Keek-Beyer) 
displacement  solutions,  and  the  difference  between  the  invis- 
cid and  viscous  (parametric  differentiation  -  P.D.)  displace- 
ment solutions,  respectively,  for  e  =  .01,  and  ij>  =  .01  (T  = 
1.27).   The  two  plots  of  the  calculated  difference  are  the  ba- 
sic results  used  for  evaluation  of  the  parametric  differentia- 
tion solution.   The  differences  are  calculated  according  to 
Equation  4.16,  and  one  would  expect  the  calculated  differences 
to  be  in  phase  with  the  displacement  solutions;  since  viscos- 
ity should  reduce  the  absolute  amplitudes,  the  difference 
should  be  maximum  when  the  Fubini  solution  has  a  maximum  and 
minimum  when  the  Fubini  solution  has  a  minimum.   The  zeroes 
will  not  necessarily  coincide  since  the  effect  of  the  non-lin- 
earity is  to  lengthen  the  compressive  portion  of  the  waveform 
while  reducing  its  magnitude  and  to  shorten  the  rarefactive 
portion  while  increasing  its  absolute  magnitude;  viscosity  will 
tend  to  keep  the  waveform  sinusoidal,  thus  maintaining  the 
length  of  the  two  portions  of  the  waveform.   Thus  one  could 
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hardly  expect  the  zeroes  to  coincide. 

Since  T  =  1.27,  one  can  infer  that  the  Keek-Beyer  solu- 
tion should  be  accurate,  and  also  that  some  harmonic  content 
should  be  evident  in  the  particle  velocity  solution,  assuming 
T  =  1  is  an  accurate  representation  of  the  incipience  of  im- 
portant non-linear  effects.   In  analyzing  Figures  5c  and  5d , 
one  is  immediately  struck  by  the  fact  that  although  the  calcu- 
lated differences  seem  to  be  of  the  same  form,  there  seems  to 
be  some  sort  of  ramp  type  behavior  entering  into  the  P.D.  dis- 
placement solution.   If  one  were  to  try  to  place  a  physical 
interpretation  on  this  result,  it  would  be  that  the  ability  of 
the  medium  to  restore  itself  elastically  has  been  exceeded  due 
to  the  large  displacement  amplitudes  and  that  therefore  a  per- 
manent displacement  has  been  induced,  just  as  if  a  spring  had 
been  stretched  beyond  its  elastic  limit  inducing  a  permanent 
strain.   Due  to  the  fluid  nature  of  the  medium  and  because  no 
analytic  solution  predicts  such  behavior,  such  a  physical  in- 
terpretation may  not  have  much  meaning.   An  attempt  will  be 
made  to  offer  a  possible  explanation  for  this  result  in  a  la- 
ter paragraph,  but  this  remains  an  unsolved  problem  in  this 
work.   The  average  numerical  loss  in  the  positive  portion  of 
the  waveform  is  clearly  greater  than  that  of  the  negative  por- 
tion, rather  than  being  equal  as  one  would  expect. 

The  velocity  solution,  on  the  other  hand,  seems  to  be 
quite  normal.   One  can  observe  that  though  non-linear  effects 
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are  producing  some  steepening  of  the  profile,  viscous  dissipa- 
tion does  in  fact  reduce  the  amplitudes  properly;  the  ampli- 
tude of  the  last  peak  is  3.1  dB  down  vice  3.3  dB  down  by  the 
Keek-Beyer  solution,  a  difference  of  about  6%.   The  last  point 
in  the  velocity  profile  is  numerically  inaccurate  because  of 
the  lack  of  sufficient  points  at  the  tip  of  the  solution  mesh. 
The  small  "bump"  in  the  velocity  profile  is  not  easily  ex- 
plained, but  probably  arises  from  numerical  inaccuracies.   It 
seems  remarkable  that  with  little  direct  effort  aimed  at  re- 
ducing computational  inaccuracies  in  the  numerical  scheme,  the 
difference  in  the  solutions  is  but  6%,  remembering  all  the 
while  that  this  difference  may  not  represent  an  error  at  all, 
but  rather  a  more  accurate  solution  produced  by  a  method  which 
does  not  resort  to  common  approximations.   One  may  wonder  at 
this  point  why  an  equation  in  particle  displacement  vice  one 
in  particle  velocity  was  treated,  the  answer  being  that  the 
equation  in  particle  velocity  would  still  have  required  one  to 
know  the  displacement  solution  in  order  to  arrive  at  the  vari- 
able coefficients  which  would  have  appeared  in  the  parameter 
space  equation,  and  thus  would  probably  have  been  no  essential 
simplification.   Another  aspect  of  the  P.D.  solution  which  is 
worthy  of  note  is  the  steepness  of  the  velocity  profile.   The 
Keek-Beyer  solution  for  this  value  of  Y    looks  essentially  like 
a  decaying  sinusoid,  while  the  P.D.  solution  appears  to  indi- 
cate that  the  non-linearity  remains  fairly  important,  result- 
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ing  in  some  steepening  of  the  profile.   This  suggests  that 
perhaps  r  =  1  may  not  be  a  sufficiently  strong  indicator  of 
the  incipience  of  non-linear  effects. 

Results  analogous  to  those  of  Figure  5  have  been  obtained 
for  different  values  of  the  problem  parameters,  but  with  the 
number  of  points  utilized  in  the  calculation  maintained  at 
N  =  .01. Figures  6  and  7  are  the  plotted  results  for  e  =  .01 
and  t\>   =    .001  (T  =  12.7),  and  for  e  =  .01  and  \p   =    .0000463  (V 
=   275)  respectively.   Figures  8,  9,  and  10  are  the  results  for 
e  =  .02  and  \\>   =    .01  (T  =  2.55),  ty   =  .001  (T    =    25.5),  and  if;  = 
.0000463  (T  =  550)  respectively.   In  addition,  calculations 
have  been  conducted  for  e  =  .0075  and  e  =  .005,  but  are  not 
shown  here.   These  results  are  not  much  different  in  form  from 
those  discussed  previously,  but  they  uncover  some  problems  with 
the  solution  technique  as  well  as  help  demonstrate  that  the 
method  is  in  fact  a  potentially  useful  one  for  investigation 
of  non-linear  wave  propagation  problems,  since  it  successfully 
solves  the  equation  of  motion  in  a  predictable  and  consistent 
manner  for  various  values  of  r.   The  range  of  the  3£  product 
had  to  be  restricted  as  it  was  because  of  the  need  to  minimize 
computer  costs  and  thus  the  number  of  points  at  which  the  dif- 
ference equation  was  to  be  solved. 

Shifting  attention  to  Figure  6,  first  recall  that  with  T 
=  12.7  the  Keek-Beyer  solution  is  probably  only  marginally  ac- 
curate.  The  displacement  solution  again  exhibits  the  ramp  type 
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behavior,  although  it  is  only  readily  clear  when  analyzing 
Figure  6d ,  which  shows  a  greater  average  loss  in  the  positive 
portion  of  the  waveform  than  in  the  negative  portion.   Figure 
6c  shows  that  the  Keek-Beyer  solution  is  in  fact  no  longer  ac- 
curate near  the  discontinuity  distance.   The  velocity  waveform 
appears  to  be  consistent  with  expectations  with  only  slight 
attenuation  because  of  the  relative  smallness  of  the  attenua- 
tion coefficient  over  the  spatial  length  of  the  calculations. 
However,  small  signal  attenuation,  i.e.,  that  attenuation  which 
would  occur  if  no  harmonics  were  non-linearly  produced,  would 
predict  a  local  Mach  number  of  .0092  at  the  discontinuity  while 
in  this  case  a  value  of  .0096  was  obtained.   This  result  is 
not  consistent;  one  would  expect  greater  attenuation  in  the 
non-linear  case  because  the  generation  of  harmonics  will  re- 
duce the  magnitude  of  the  primary  in  addition  to  small  signal 
attentuation,  while  the  harmonics  themselves  will  be  attenu- 
ated even  more  rapidly  because  of  greater  attenuation  coeffi- 
cients.  Unfortunately,  in  this  case  there  is  no  easy  explana- 
tion of  the  result;  however,  it  did  prompt  some  thoughts  about 
what  the  Fubini  solution  looks  like  and  what  the  changing 
shape  of  the  waveform  should  imply  with  respect  to  the  instan- 
taneous values  of  particle  velocity  near  the  discontinuity 
distance. 

From  the  point  of  view  of  energy  conservation,  one  would 
expect  that  in  a  fluid  without  viscous  dissipation,  the  time 
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or  space  average  energy  density  in  a  plane  wave  train  propagating 
isentropically  will  remain  constant.   The  implication  of  this 
is  that  as  the  finite  amplitude  waveform  steepens  and  eventu- 
ally assumes  a  sawtooth  shape,  the  peak  values  of  the  particle 
velocity  must  increase  in  order  to  maintain  constant  average 
energy  density  in  the  wave  train.   For  a  sawtooth  wave  to  have 
the  same  energy  density  as  a  sinusoid  of  unit  amplitude  and 
identical  period,  its  amplitude  must  be  approximately  1.22. 
Realizing  that  a  fully  formed  sawtooth  develops  at  a  distance 
from  the  boundary  of  about  3 . 5X  [3],  it  would  still  be  expected 
that  the  peak  value  of  the  particle  velocity  at  the  discontin- 
uity distance  would  be  greater  than  that  at  the  boundary.   The 
base  solution  (Fubini)  calculated  in  this  study  is  accurate  to 
five  places  in  displacement  and  three  in  velocity,  and  the 
peak  particle  velocity  at  the  discontinuity  distance  is  ident- 
ical to  that  at  the  boundary,  implying  that  there  in  fact  has 
been  an  energy  loss  in  an  inviscid  solution,  and  it  is  certain- 
ly not  clear  where  this  energy  has  gone.   In  any  case,  for  a 
T  of  12.7,  non-linear  effects  are  definitely  important,  and  it 
is  entirely  conceivable  that  the  P.D.  solution  has  somehow 
recognized  an  error  in  the  base  solution  and  has  yielded  a 
quadrature  value  which  correctly  accounts  for  the  steepening 
of  the  profile  and  thus  predicts  a  particle  velocity  which  is 
greater  than  that  which  would  be  obtained  with  a  sinusoid  and 
small  signal  attenuation.   This  explanation  is  of  course  sup- 
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position  since  there  are  too  few  points  calculated  to  obtain 
a  good  space  average  energy  density  to  verify  an  energy  loss 
which  is  greater  than  that  predicted  by  small  signal  attenua- 
tion.  It  indicates  an  unsolved  problem  and  at  the  same  time 
suggests  that  parametric  differentiation  predicts  the  proper 
solution  to  the  problem.   Blackstock  [9]  suggests  the  mechan- 
ism for  dissipation  in  an  inviscid  fluid  after  a  shock  has 
formed;  here,  there  is  no  reason  for  such  dissipation  since 
there  are  no  shocks  in  the  first  zone  of  non-linear  propaga- 
tion. 

Figure  7  shows  much  the  same  results  as  those  of  Figures 
5  and  6.   In  this  case,  r  =  275,  and  the  Keek-Beyer  solution 
is  clearly  not  converging  to  the  proper  result  as  evidenced  by 
Figure  7c.   Again  the  ramp  behavior  is  clear  in  Figure  7d, 
though  the  magnitudes  of  the  difference  values  are  understand- 
ably smaller.   The  peak  value  of  the  calculated  Mach  number  at 
the  discontinuity  distance  is  now  .0106  which  is  in  fact  great- 
er than  the  peak  at  the  boundary.   The  differences  between  the 
inviscid  and  viscous  displacement  solutions  are  now  so  small 
that  on  the  scale  of  these  plots,  the  viscous  solution  appears 
identical  to  the  inviscid  solution.   This  is  not  remarkable, 
but  it  is  noteworthy  that  the  plotted  differences,  though 
small,  produce  a  smooth  curve  rather  than  one  that  jumps  from 
point  to  point.   Figures  8,  9,  and  10,  which  represent  the 
calculations  for  a  source  Mach  number  of  .02,  show  results 
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which  are  similar  to  those  shown  in  Figures  5,  6,  and  7,  and 
do  not  need  to  be  analyzed  with  as  much  detail.   However,  in 
Figure  8,  with  r  =  2.55,  the  peak  Mach  value  near  the  discon- 
tinuity distance  is  .0147  vice  .0137  according  to  the  Keek- 
Beyer  solution  for  a  difference  of  7%,  which  is  consistent 
with  the  6%  difference  obtained  in  the  case  e  =  .01,  T   =    1.27. 
It  is  also  consistent  with  the  idea  that  peak  Mach  numbers 
should  increase  in  the  non-linear  case  since  small  signal  at- 
tenuation would  predict  a  Mach  number  of  .014  at  the  same 
point  for  a  purely  sinusoidal  wave.   The  ramp  behavior  is  evi- 
dent in  each  of  the  Figures  8d,  9d ,  and  lOd. Finally,  these  results 
also    provide  a  counter  example  for  the  idea  that  the  peak 
Mach  number  should  increase  near  the  discontinuity  distance  in 
the  case  of  strong  non-linearity,  i.e.  when  V    is  greater  than 
10.   The  particle  velocity  plots  of  Figures  9b  and  10b  indi- 
cate peak  Mach  values  which  remain  less  than  the  peak  at  the 
boundary;  in  Figure  9b  the  value  is  .0193  and  in  Figure  10b  it 
is  .01998.   Small  signal  attenuation  would  predict  a  value  of 
.0192  for  the  result  of  Figure  9b.   Energy  considerations  dic- 
tate an  increase  in  the  peak  Mach  values  and  in  this  case  it 
does  not  occur.   However,  the  spatial  step  size  used  to  arrive 
at  the  plots  of  Figures  6  and  7  was  .04,  while  in  the  calcula- 
tions used  to  obtain  the  plots  of  Figures  9  and  10,  the  spa- 
tial step  size  was  .02.   Recalling  the  finite  difference  ap- 
proximations utilized,  this  difference  in  step  size  may  be  the 
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explanation  for  the  difference  in  results.   The  truncation  er- 
ror for  the  second-order  difference  approximations  is  of  OuAx)2.] 
Therefore,  it  is  certainly  conceivable  that  the  step  size  of 
.04  is  not  accurate  enough  to  yield  good  results  for  particle 
velocity.   At  the  same  time  the  considerations  of  energy  can- 
not be  ignored.   It  is  likely  that  the  increase  in  the  peak 
Mach  number  in  the  results  of  Figures  6  and  7  is  a  result  of 
numerical  inaccuracy,  but  this  is  certainly  not  entirely  clear. 
The  other  calculations  conducted  with  Mach  numbers  of  .007  5 
and  .005  produced  results  which  were  consistent  with  the  other 
results  obtained,  but  they  are  not  presented  here  because  the 
spatial  step  size  was  sufficiently  large  that  the  results  are 
somewhat  questionable  because  of  truncation  error. 

There  remains  a  question  with  respect  to  the  curious  ramp 
behavior  in  the  particle  displacement  solutions.   In  examining 
the  plots  of  the  difference  between  the  inviscid  and  viscous 
solutions,  this  behavior  is  always  clear,  but  would  not  neces- 
sarily be  clear  at  all  in  analyzing  the  particle  displacement 
results,  particularly  as  ip  becomes  smaller  and  smaller.   i|>  is 
a  direct  measure  of  the  viscosity  of  the  assumed  fluid  rather 
than  a  relative  measure  of  viscosity  such  as  T.      Recalling  that 
a  i|j  of  .01  would  imply  a  highly  viscous  liquid  or  a  very  high 
frequency  wave  and  that  the  approximation  upon  which  previous 
authors1  results  are  based  is  a/K  <<  1,  which  implies  weak 
viscosity,  the  reason  for  the  ramp  may  become  clear.   Previous 


-66- 


analyses,  in  making  this  approximation,  may  have  neglected  a 
rather  interesting  physical  interpretation  of  the  non-linear 
loss  mechanism.   For  small  \\> ,    the  values  of  the  calculated 
differences  are  so  small  compared  to  the  values  of  particle 
displacement  that  this  effect  would  not  be  noted   when  solv- 
ing for  displacement  or  velocity  without  consideration  of  the 
differences  between  the  viscous  and  inviscid  solutions.   In 
fact,  in  making  the  assumption  a/K  <<  1,  the  values  of  ty   are 

apparently  restricted  to  relatively  small  values;  a/K  =  ^ — j- 

Tvb      ,  ,  vb    ,  ,    -     •  .  •    t     .  t_  j_ 

=  - — v  ;  and  iL  =  - — r-  ,  therefore  it  is  clear  that  previous 
LA  C  A 

o  o 

analyses  have  not  considered  values  of  ip  on  the  order  of  those 
considered  here.   For  large  ip  and  thus  highly  viscous  propaga- 
tion, previous  results  are  not  valid,  and  therefore  if  the 
ramp  behavior  has  a  physical  explanation,  it  could  not  have 
been  noted  before. 

In  the  parameter  space  Equations  4.6  or  4.9,  there  are 
four  terms  which  involve  the  dependent  variable/   g   ,  g   , 

XX      X- L- 

g    ,  and  g  t       which  result  from  similar  terms  in  the  displace- 

XXX-  X 

ment  space  Equation  3.9.   g    is  a  term  which  relates  to  the 
elastic  properties  of  the  medium,  g    to  the  inertial  proper- 
ties, and  g    to  the  viscous  properties.   The  g   term  implies 

XX  X-  X 

that  there  exists  some  modification  of  the  elastic  properties 
of  the  medium.   Recalling  the  discussion  in  the  previous  sec- 
tion concerning  the  effect  of  this  term  in  the  portions  of  the 
waveform  where  the  velocity  gradient  u   assumes  large  values, 
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the  source  of  the  ramp  behavior  may  become  evident.   When  E, 
is  large,  the  g   term  will  dominate  the  equation.   In  the  in- 

viscid  solution,  E        becomes  infinite  at  the  first  shock. 

xx 

Since  the  inviscid  solution  is  in  fact  the  starting  solution 
and  the  quadrature  from  parameter  space  to  displacement  space 
is  accomplished  with  an  integration  with  respect  to  \|j  ,  it  is 
obviously  not  unreasonable  to  expect  E,        to  assume  large  val- 
ues  when  \p    is  small;  also,  one  would  expect  E,        to  assume  nu- 
merical  values  which  are  more  accurate  for  the  current  \\)    step 
of  the  calculation  if  the  spatial  step  size  was  small  in  the 
regions  of  large  E,      .   In  any  case,  if  the  g   term  dominates 
the  equation  in  a  specific  region,  the  g  solution  in  that  re- 
gion will  vary  linearly  with  x,  hence  the  E,    solution  will  vary 
in  a  similar  manner,  and  hence  the  ramp. 

In  the  absence  of  viscosity,  it  is  clear  that  the  fluid 
particles  oscillate  about  an  equilibrium  position;  with  vis- 
cosity, realizing  that  the  g   term  is  one  which,  in  essence, 
alters  the  elastic  constant  of  the  system,  it  is  postulated 
that  the  following  phenomenom  occurs.   Because  of  the  large 
velocity  amplitudes,  the  phase  velocity  increases  in  the  com- 
pressive portion  of  the  wave  and  decreases  in  the  rarefactive 
portion;  this  is  well  known.   These  large  particle  velocities 
produce  correspondingly  large  particle  displacements,  but  now 
the  fluid  particles,  instead  of  being  purely  elastic,  do  not 
have  all  of  their  kinetic  energy  at  zero  displacement  converted 


-68- 


to  potential  energy  at  full  displacement.   Some  is  lost  to 
small  signal  attenuation   of  both  fundamental  and  harmonics, 
and  the  rest  produces  an  irreversible  displacement  of  the 
fluid  particles;  that  is,  work  is  expended  on  the  fluid  system 
to  produce  a  displacement  of  the  equilibrium  positions  of  the 
fluid  particles.   Unfortunately,  this  effect  can  only  be  pos- 
tulated here  because  a  precise  energy  balance  has  not  been 
performed.   In  order  to  calculate  the  specific  energy  lost  in 
producing  a  given  displacement  at  a  given  point,  the  net  ac- 
celeration on  the  fluid  particles  which  produced  this  displace- 
ment would  be  required,  and  this  calculation  has  not  yet  been 
performed. 

The  reason  for  this  apparent  displacement  of  the  equili- 
brium positions  of  the  fluid  particles  may  be  explained  in  a 
more  rhetorical  manner  as  an  interaction  between  non-linear 
and  viscous  effects.   The  non-linear  effect  is  to  alter  the 
elastic  constant  thus  changing  the  speed  of  propagation  of  the 
wave,  resulting  in  the  steepening  of  the  waveform.   Now,  be- 
cause the  oncoming  compressive  portion  of  a  given  cycle  of  the 
wave  is  travelling  faster  than  the  preceding  rarefactive  por- 
tion and  because  viscosity  is  acting  to  slow  the  velocity  of 
the  fluid  particles  in  the  rarefactive  portion  even  more,  the 
particles,  after  passage  of  the  given  cycle,  do  not  return  to 
their  former  equilibrium  position,  but  are  instead  minutely 
displaced.   The  oncoming  compressive  portion  of  the  wave  is 
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' filling1  the  preceding  raref active  portion  and  consequently 
the  fluid  particles  in  the  rarefactive  portion  do  not  have  the 
opportunity  to  recover  their  initial  equilibrium  positions. 
Such  an  effect  would  obviously  tend  to  create  a  vacuum  behind 
the  wave  train,  but  this  can  be  justified  by  requiring  that 
fluid  flow  in  from  outside  the  region  excited  by  the  source. 

Finally,  there  are  some  matters  concerning  the  numerical 
method  which  require  attention.   The  stability  criterion  de- 
rived in  the  preceding  chapter  was  arrived  at  in  a  rather  un- 
conventional manner.   Some  numerical  experimentation  was  per- 
formed to  determine  just  how  sensitive  to  the  ratio  of  step 
sizes  the  solution  was.   The  exact  stability  criterion  was 


Ax  <  At/(1  +  |^)2e  (6.1) 


which  implies  for  9£/8x  positive,  e  =  .01,  and  3=4  that  Ax 
<_   At/1.08  near  the  first  shock.   This  should  be  an  absolute 
limit  on  the  spatial  step  size  since  8£/9x  will  have  its  maxi- 
mum value  at  this  point.   The  equations  were  solved  for  Ax=.04 
and  At=.04,  and  the  numerical  solution  became  highly  unstable 
beyond  the  second  cycle  of  the  waveform.   Then  At  was  set  equal 
to  .045  which  would  be  just  slightly  greater  than  the  At  re- 
quired by  the  stability  criterion  for  Ax=.04,  and  the  solution 
was  stable  throughout.   All  subsequent  calculations  were  con- 
ducted with  a  convenient  At  selected  so  that  it  was  large  enough 
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for  the  given  Ax  to  satisfy  the  stability  criterion. 

Another  matter  of  concern  in  the  calculations  is  the  ef- 
fect of  the  step  sizes,  in  x  and  t  and  in  ip  as  well,  on  the 
accuracy  of  the  solution.   In  this  analysis,  the  base  solution 
was  calculated  accurate  to  five  significant  figures,  obviously 
limiting  the  accuracy  of  the  calculated  solution  to  the  same 
extent.   The  base  solution  required  51  harmonics  at  the  dis- 
continuity distance  for  this  accuracy.   The  spatial  step  sizes 
used  for  the  various  solutions  were  Ax  =  .08,  .06,  .04,  and 
.02.   Since  the  truncation  errors  in  the  second  order  differ- 
ences are  of  OuAx)  2J,  it  would  be  reasonable  to  assume  that  the 
calculated  values  of  g  would  be  accurate  to  two  or  three  places 
depending  on  the  step  size  used.   However,  to  return  to  dis- 
placement space,  the  values  of  g  are  multiplied  by  a  difference 
between  two  values  of  \p ,    and  hence  for  a  maximum  ty   of  .01,  if 
g  is  accurate  to  three  significant  figures,  the  differences 
calculated  in  the  quadrature  will  leave  the  displacement  solu- 
tion accurate  to  five  significant  figures,  neglecting  errors 
introduced  by  the  discretization  of  the  solutions  with  respect 
to  ip.   This  then  would  be  consistent  with  a  base  solution  cal- 
culated to  four  significant  figures,  for  Ax=.02  and  marginally 
consistent  for  Ax=.04.   This  is  the  reason  why  the  results  for 
e  =  .005  and  e  =  .0075  were  not  presented;  their  accuracy  could 
only  be  considered  marginal  at  best  for  the  maximum  ip  value. 
These  were  the  solutions  which  employed  step  sizes  of  .06  and 
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.08.   Roundoff  error  was  not  considered  to  be  of  significant 
concern  in  this  analysis  with  machine  accuracy  to  seven  signi- 
ficant figures.   The  solutions  also  did  not  seem  to  be  sensi- 
tive to  the  step  size  in  time  as  long  as  the  time  step  was 
large  enough  to  satisfy  the  stability  criterion  at  all  points. 
Finally,  the  step  size  in  \p    seemed  to  matter  very  little  in 
the  form  of  the  results  if  not  in  the  exact  numerical  results. 
It  is  obvious  that  in  integrating  across  orders  of  magnitude, 
the  integration  must  be  handled  carefully,  but  experience  dem- 
onstrated that  values  of  g  at  a  given  point  changed  little  as 
ty   was  varied,  and  thus  the  error  introduced  by  utilizing  a 
relatively  large  step  size  in  ip  should  be  small.   Consequently, 
the  following  scheme  was  adopted: 

1)  the  first  ty   value  was  the  limiting  value  ip 

2)  the  second  ip  value  was  that  which  when  used  in 
the  quadrature  yielded  difference  values  which 
were  just  within  machine  roundoff  accuracy 

3)  each  subsequent  step  in  the  order  of  magnitude 
of  \p   was  broken  into  an  algebraically  increasing 
number  of  steps  on  a  logarithmic  basis. 


—  f> 
For  example,  with  tJj   =  0,  if>,  =  10   ,  there  would  be  two  steps 

-6  -5  -5 

between  \p   =  10    and  ij;  =  10   ,  with  \\>?   =    .316  x  10    and  iK 

-5  -5  -4 

=  10   ,  three  steps  between  ip  =  X 0    and  ty   -   10   ,  4>.  =  .214 

-4  -4  -4 

x  10  ,    ty-    =    .463  x  10   ,  and  ib,  =  10   ,  etc.   There  is  no 
~>  b 
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clear  method  to  determine  the  proper  step  size  in  \p ,    but  this 
one  produced  consistent  results.   In  general,  it  would  be  in- 
appropriate if  g  varied  in  a  highly  non-linear  manner  with  \p . 
The  final  comments  concerning  the  numerical  results  are 
made  with  respect  to  the  formulation  of  the  problem.   Equations 
4.6  and  4.9  were  the  parameter  space  equations  developed  from 
the  equation  treated  by  Keck  and  Beyer  and  from  the  parent  e- 
quation  without  approximation,  respectively.   Early  calcula- 
tions were  performed  with  Equation  4.6  and  later  ones  with  E- 
quation  4.9.   For  e  =  . 01  and  ty   =  .01,  the  maximum  percentage 
difference  between  the  calculated  solutions  by  Equations  4.6 
and  4.9  was  approximately  4%,  demonstrating  that  the  equation 
solved  in  previous  work  does  in  fact  produce  very  little  error 
for  Mach  numbers  up  to  .01.   Finally,  the  two  forms  of  the 
boundary  values  for  g  appeared  to  make  little  difference  in 
the  final  result;  producing  a  maximum  percentage  difference  in 
the  calculated  values  of  particle  displacement  of  approximate- 
ly 1%  between  the  case  utilizing  the  'natural'  boundary  condi- 
tions for  g  and  the  case  in  which  artificial  dependence  on  ij; 
was  introduced  at  the  boundary. 
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VII.  CONCLUSIONS  AND  RECOMMENDATIONS 

A  theoretical  study  of  non-linear  plane  acoustic  wave 
propagation  in  a  viscous  fluid  has  been  conducted  utilizing  a 
little  known  mathematical  technique,  parametric  differentiation, 
in  order  to  demonstrate  its  utility  when  applied  to  acoustic 
propagation  problems,  and  to  obtain  whatever  physical  insight 
into  the  problem  the  results  may  present.   It  has  been  demon- 
strated that  the  method  can  in  fact  yield  results  which  are 
consistent  with  other  analytic  solutions  to  the  problem  for 
certain  arbitrarily  selected  values  of  the  viscous  and  non- 
linear parameters.   No  restriction  was  required  in  the  selec- 
tion of  these  values  by  the  mathematical  nature  of  the  method, 
and  there  is  no  reason  to  believe  that  the  method  could  not  be 
applied  for  any  values  of  the  parameters.   Restrictions  which 
may  exist  in  the  application  of  the  method  are  computational 
ones,  specifically:   1)  the  results  can  only  be  as  accurate  as 
the  solution  to  the  problem  for  ip  =  \p    ;  2)  accuracy  may  be 
further  limited  by  the  finite  difference  approximations  util- 
ized; 3)  accuracy  may  be  limited  by  the  ability  of  the  computer 
to  solve  the  finite  difference  scheme  for  a  large  number  of 
points  with  a  given  time  and  memory  budget;  for  instance,  dou- 
ble precision  will  require  twice  the  computer  memory;  and  4) 
the  hyperbolic  nature  of  the  problem  requires  calculations  to 
be  conducted  throughout  the  domain  of  dependence  of  a  given 
point  to  arrive  at  a  solution  at  that  point.   These  restric- 
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tions  may  render  the  method  unusable  for  calculations  involv- 
ing large  propagation  distances,  but  they  do  not  detract  from 
its  ability  to  provide  physical  insight.   Further  research  may 
provide  a  reasonable  technique  for  application  of  the  method 
when  large  distances  are  involved.   The  solution  of  the  prob- 
lem by  parametric  differentiation  has  allowed  the  elimination 
of  the  restriction  a/k  <<  1  heretofore  assumed  in  typical  pro- 
pagation studies.   Elimination  of  this  restriction  has  produced 
results  which  are  valid  for  highly  viscous  fluids  as  well  as 
demonstrated  that  the  results  for  such  fluids  are  essentially 
similar  to  those  for  less  viscous  fluids.   It  has  allowed  the 
formulation  of  a  postulated  behavior  for  all  viscous  fluids, 
which  may  explain  physically  why  there  are  non-linearly  in- 
duced energy  losses  above  and  beyond  those  associated  with 
viscous  attenuation.   Specifically  the  solutions  produced  by 
this  method  predict  that  some  of  the  energy  contained  in  a 
finite  amplitude  wave  in  a  viscous  medium  is  expended  in  pro- 
ducing a  displacement  of  the  equilibrium  positions  of  the  fluid 
particles.   These  results  have  also  been  obtained  cheaply;  for 
101  spatial  points,  the  computer  time  required  to  produce  a 
result  for  single  values  of  ty   and  the  $e  product  is  approxi- 
mately .15  minutes  with  400  kilobytes  of  computer  memory. 

This  is  not  to  suggest  that  this  work  proves  beyond  a 
doubt  that  the  method  is  a  useful  one.   Consideration  of  prop- 
agation once  a  shock  has  actually  formed  has  not  been  attempted 
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If  the  shocks  were  in  fact  discontinuities  as  predicted  by  the 
inviscid  theory,  the  technique  probably  would  not  be  useful 
beyond  the  first  zone  of  propagation.   However,  since  the  wave- 
form is  in  fact  continuous  through  a  shock  of  finite  thickness, 
it  is  reasonable  to  expect  that  with  denser  sampling  of  the 
waveform  through  the  shocks,  good  numerical  results  could  be 
obtained.   The  next  step  in  a  continuing  study  of  the  solution 
of  the  problem  through  parametric  differentiation  would  be  to 
select  Blackstock's  weak  shock  solution  [3],  which  connects 
the  Fubini  solution  [1]  with  the  Fay  solution  [2]  in  its  in- 
viscid  limit  to  yield  a  solution  valid  to  the  end  of  the  sec- 
ond zone  of  propagation,  as  the  base  solution  for  i|j   =0,  and 
to  attempt  to  extend  the  solution  into  the  region  of  periodic 
shocks.   With  this  result  a  precise  energy  balance  calculation 
would  demonstrate  whether  or  not  energy  is  lost  in  producing  a 
permanent  displacement  of  the  fluid  particles. 

More  fundamentally,  it  would  be  useful  if  a  mathematical 
analysis  of  the  technique  were  to  be  conducted  in  order  to  de- 
termine why  it  apparently  correctly  predicts  solutions  for  many 
different  classes  of  problems,  and  what  if  any  are  the  theor- 
etical limitations  to  its  application.   Beyond  this,  there  are 
more  interesting  problems  in  non-linear  wave  propagation  to 
which  the  method  could  be  applied.   The  first,  and  one  which 
can  be  readily  adapted  to  the  computer  program  already  written, 
would  be  an  investigation  of  the  solution  of  the  Burgers  equa- 
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tion  as  Y    is  varied.   Having  noted  previously  that  this  solu- 
tion is  difficult  to  calculate  as  r  becomes  much  greater  than 
x/X  or  for  r  >  50  because  of  slow  convergence,  it  is  suggested 
that  the  analytic  solution  be  calculated  for  r  =  1.   Then  with 
the  Burgers  equation  in  the  form  solved  by  Blackstock  [5] , 

with  V  =  u/U  , 
o 


V   -  VV   =  r  1V   ,  V(o,  y)  =  sincoy  (7.1) 

O        y  yy '       '  J  J  ' 


differentiate  with  respect  to  r  to  obtain, 


g   -  Vg   -  V  g  =  r~1g    -  T    2V   ,  g(o,  y)  =  0        (7.2) 

^g    ^y    y^  yy      yy 


where  a  is  a  spacelike,  and  y  a  timelike  variable.   This  equa- 
tion would  be  perfectly  suited  to  solution  on  the  mesh  used  in 
this  analysis  because  there  would  be  no  value  at  the  i+l,j-l 
point,  and  the  scheme  would  be  completely  explicit.   The  solu- 
tion for  T  =  1  could  then  be  extended  to  arbitrary  values  of  r with- 
out  significant  computational  problems.   Because  of  the  space 
scale  of  Equation  7.1,  this  technique  will  allow  easy  calcula- 
tion of  the  solution  in  both  the  first  and  second  zones  of 
propagation. 

Finally,  there  are  the  two  problems  of  finite  amplitude 
wave  propagation  which  are  yet  to  be  solved  analytically,  those 
of  propagation  of  cylindrical  and  spherical  waves  in  viscous 
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fluids;  parametric  differentiation  may  provide  an  invaluable 
tool  in  solving  these  problems.   The  fundamental  equations  of 
motion  do  not  lend  themselves  to  analytic  solution,  there  be- 
ing no  simple  transform  such  as  that  used  by  Blackstock  in  the 
plane  wave  problem  to  reduce  the  equation  of  motion  to  a  Bur- 
gers equation  [18].   Indeed,  the  inviscid  forms  of  the  equa- 
tion have  no  exact  analytic  solutions  which  have  been  published 
Blackstock  [9]  has  obtained  approximate  solutions  to  the  cylin- 
drical and  spherical  wave  problems  for  an  inviscid  fluid,  one- 
dimensional  wave  propagation,  and  r  >>  1,  which  are  very  sim- 
ilar in  form  to  the  Fubini  solution;  these  results  were  ob- 
tained through  manipulations  of  the  equation  of  motion  which 
yielded  an  inviscid  Burgers  equation.   Such  solutions  could  be 
used  as  one  base  solution  for  a  solution  to  the  viscous  prob- 
lem by  parametric  differentiation,  realizing  that  they  are  on- 
ly approximate.   Other  authors  [12] [19]  have  produced  particu- 
lar solutions  for  spherical  and  cylindrical  wave  propagation, 
but  with  an  accurate  base  solution,  parametric  differentiation 
could  be  used  to  generate  many  solutions,  for  a  fixed  3e  pro- 
duct, with  relative  economy.   Consequently,  the  apparent  course 
of  study  to  pursue  for  non-linear  cylindrical  and  spherical 
wave  propagation  is  the  development  of  improved  base  solutions 
for  \\>      =0  and  the  application  of  parametric  differentiation 
to  these  solutions,  producing  a  large  number  of  solutions  which 
may  then  allow  for  some  generalization  of  the  results.   Such  a 
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research  program  could  lead  to  a  significant  gain  in  the  under- 
standing of  the  physical  mechanisms  which  are  collectively  des- 
cribed as  non-linear  attenuation  in  viscous  fluids. 
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Figure  2 The  Solution  Procedure 
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Figure  2 The  Solution  Procedure 
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Figure    8d.  Difference   Between   Fubini   and   P.D 

Displacement   Profiles,    3e    =    .08,    r    =    2.55 
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